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Preface 


The study of crystal and its growth has a long tradition. The art of crystal growth is 
now well developed to manufacture various high quality crystals for electronics, but 
its scientific understanding is still developing. The modern technologies and obser¬ 
vation methods of the atomic scale make possible and also require the microscopic 
understanding of the growth mechanisms to control the quality of the product. 

Since the crystal growth is a typical example of a system far from equilibrium 
which involves the first-order phase transition, its comprehension involves many dis¬ 
ciplines. In these two decades novel theoretical contributions are given from the field 
of statistical physics: Kosterlitz-Thouless phase transition of the surface roughness 
in equilibrium surface structure, and the microscopic solvability criterion for the se¬ 
lection of the dendritic crystal growth. 

This note tries to give the concise and precise overview on the fundamental scien¬ 
tific aspects of crystal growth from the simplest phenomenological argument to the 
latest discoveries. Technology and art of crystal growth are omitted in this book. 
Important coupling of the hydrodynamics to the crystal growth are beyond the ca¬ 
pability of the author and arc not touched upon. 

The note is based on my introductory lectures on the fundamental and physico- 
mathematical aspects of the crystal growth given at a winter school of Japanese 
Association of Crystal Growth and in graduate courses in many universities in Japan; 
Keio, Ochanomizu, Hokkaido, Kyoto and so on. 

I am especially grateful to H. Muller-Krumbhaax. From discussions and a long- 
lasting collaboration with him, I learned most of the fundamental theories of crystal 
growth described in this book. The collaborations and enlightening discussions with 
C. Misbah, E. Brener, D. Temkin and M. Uwaha are also deeply acknowledged. 
To M. Uwaha I am thankful for his reading through and giving comments on the 
manuscript. I am also indebted to T.Ohta, K.Wada, M.Kitamura, who have invited 
me to their universities and given chances to teach the courses which eventually 
resulted in this book. Lastly I dedicate my heartful appreciation to late Prof. R. 
Kubo, who introduced me to the fascinating world of statistical physics. 


Yukio Saito 
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1 Introduction 


We axe living in the age of informatics; In computers a mass of information is handled 
in high speed by integrated circuits built on a semiconductor crystal. By means of 
a laser light emitted from semiconductor crystals, information written in a compact 
disc at high density is read out. Through the optical fiber runs information as the 
optical signal with high density and in the extremity of speed. In this information 
technology, semiconductor crystals with high quality are required. There are also 
many other crystalline materials utilized in our life; steal for cars, quartz in watch, 
sugar and salt and so on. To fabricate crystals under our control we have to know 
the mechanism of crystal growth. 

There are also various crystals grown in nature; minerals, gems as diamonds, snow 
etc. They contain information about the history how and the environment in which 
they are grown: Snow is called a letter from the sky, diamond a letter from deep in 
the earth. In order to understand their messages we have to know the dynamics of 
crystal growth. 

Crystal has an ordered arrangement of atoms or molecules in microscopic scale. 
The microscopic regularity shows up in the symmetry as is observed in various diffrac¬ 
tion patterns. Also the regular arrangement of atoms brings about the symmetry in 
crystalline shape: Some crystals take simple forms as polyhedra as shown in Fig. 1.1, 
reflecting their symmetry. Some crystals have complicated forms as dendrite, as 
shown in Fig. 1.2 and in Fig. 1.3. They are complex in the sense that all the snow 
flakes, for example, look similar with six arms but none of them are completely the 
same. We want to know how the crystal shapes are determined. The difference in 
crystal shape should be brought about by the difference in the controlling mecha¬ 
nisms of the growth dynamics. The relation between the growth mechanism and the 
resulting growth morphology is to be explored. 

There are many textbooks and monographs on crystal growth [74, 48, 120, 85]. I 
intend to give in this book a concise but precise overview on the fundamental theories 
on crystal growth from the viewpoint of statistical physics, especially on the recent 
developments in the pattern formation in the diffusion field. In Part I, the ideal growth 
formulae are derived from a thermodynamical point of view. ’’Ideal” here means that 
all the thermodynamic driving force for the phase transition is invested for the crystal 
growth. The growth formula thus gives the maximum growth velocity. In reality there 
are many hindrances against the growth. The largest effect takes place at the crystal 
surface, where the growth process takes place. In Part II, the equilibrium structure 
of the crystal surface is studied microscopically, and the surface phase transition of 
roughening is discussed. For an atomically smooth surface, surface kinetics governs 
the crystal growth and the growth laws in this situation is discussed in Part III. For 
an atomically rough surface, material transport or the heat conduction controls the 
crystal growth. In this case, macroscopic morphology of the crystal surface is strongly 
influenced, and complex pattern formation is induced. The topic will be discussed in 
Part IV. 
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1. Introduction 



Figure 1.3: Irregular dendrites of (a) MnC >2 (Courtesy by N.Osada) and (b) Au on 
Pt(lll) [41], 


/V 




Part I 

Ideal Growth Laws 
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Crystal growth is an example of a dynamical first order phase transition. A stable 
phase, crystal, grows out from a metastable phase, melt or vapor. The driving force 
for the growth is the chemical potential difference of the stable and the metastable 
phases. The simple assumption of the linear response that the growth velocity is 
proportional to the driving force gives an ideal linear growth laws. When the crystal 
is finite, or when the interface deforms, the surface tension has to be included in the 
thermodynamics of the crystal growth. Surface tension plays an important role in 
the determination of the shape of a finite sized crystal, or the deformation of the flat 
interface. 


2 Melting and Solidification: First-Order Phase 
Transition 

At low temperatures almost all the materials order in crystalline form where atoms 
are arranged regularly. (Fig.2.1a). This is the configuration with the minimum in¬ 
teraction energy E of atoms or molecules. At a high temperature atoms break the 
regular structure since it is less free and with little entropy. In an irregular arrange¬ 
ment as in liquid (Fig.2.1b) and in gas (Fig.2.1c), molecules move more freely and 
thus gain the entropy in spite of the energy cost. For a given temperature T and 
pressure p, the equilibrium state is determined by the second law of thermodynamics 


9 9 



(a) (b) (c) 


Figure 2.1: Arrangement of atoms (a) in a crystal, (b) in a liquid and (c) in a gas 
phase. 
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Part I. Ideal Growth Laws 



Figure 2.2: Equilibrium phase diagram Figure 2.3: Isobaric variation of the 
in the temperature T and the pressure Gibbs free energies of the liquid, Gl, 
p phase space. and of the crystal Gs as functions of 

temperature. 

[121] such that the Gibbs free energy G = E — TS + pV is minimized. Here V is 
the volume of the system. At low temperatures with small values of T the entropy 
S contributes less than the enthalpy H = E + pV and the crystal with low H is 
realized (Fig.2.2). At high temperatures, on the other hand, the entropy S gives the 
dominant contribution and the configuration with a large S has the low free energy 
G: the liquid or gas phase is realized. Two phases can coexist on a coexistence curve. 

We consider in the following a liquid-crystal phase transition as an example. Iso¬ 
baric variation of the Gibbs free energies of liquid and crystal, G t , and Gs respectively, 
are shown in Fig.2.3 as functions of temperature T. At a low temperature Gl lies 
higher than Gs showing that the crystal is stabler than the liquid, but on increasing 
the temperature, G L decreases faster than G s due to the large entropy Sr, of liquid 
compared to Ss of crystal. At a melting point T M (p), G r , and Gs cross with each 
other 

Gs(T M ,p) = G],(T M ,p), (2.1) 

and the liquid becomes stabler for T >T M . 

On heating the crystal under a constant pressure, its temperature first increases, 
as shown in Fig.2.4. At a melting point T M the temperature stops increasing and the 
applied heat is consumed to change the state of the matter from crystal to liquid. Since 
the absorbed heat does not appear explicitly as a temperature rise, it is called the 
latent heat L. The first law of thermodynamics says that the applied heat changes into 
the work pdV done to the environment and the increment of the internal energy d E. 
Since the pressure is kept constant, the heat changes into the enthalpy H = E + pV : 
d H = (\E + d (pV) = dE + pdV. Thus the latent heat observed at a melting point 
corresponds to the enthalpy difference of the two phases as L = H\,{Tu , p) - f/ s (T M , p), 
with Ho, and H\, being the enthalpy of the crystal and liquid phases respectively. 




3. Crystal Growth from the Melt 
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Figure 2.4: Isobaric variation of the temperature due to the heating, for example, by 
an electric current I under the voltage V. At the melting temperature T M , crystal 
starts to melt, and until the completion of melting, system absorbs the latent heat L. 

Since Gibbs free energies of both phases are equal at the melting point, G$ = Hs — 
TmSs = G l = Hi — T m S[,, the latent heat is proportional to the entropy difference 
AS = Si - S s : 

L = H l -H s = T m AS. (2.2) 

According to thermodynamics [121], entropy S is the temperature derivative of the 
Gibbs free energy as S = —(dG/dT) p . Thus the phase transition with a latent heat 
is associated with a discontinuity in the slope of the Gibbs free energy 



Ehrenfest named this type of phase transition with a discontinuity in the first deriva¬ 
tive of G as a first-order phase transition . 

3 Crystal Growth from the Melt 

Now let us cool the liquid to a temperature T below the melting point T M . The 
Gibbs free energy of a crystal G s (T,p) is lower than that of a liquid G L (T,p) as 
shown in Fig.2.3, and the true equilibrium state is a crystalline phase. However, the 
whole liquid cannot instantaneously turn into the crystal. We often experience that 
the liquid is supercooled for a long while near the melting point. Eventually, a small 
crystalline nucleus is formed in the liquid, and then it grows. The crystal grows by the 
advancement of a crystallization front. The evolution is driven by the second law of 
thermodynamics so as to minimize the Gibbs free energy at a given temperature and 
pressure [121]. The driving force of the crystal growth is the difference of the Gibbs 
free energies of the liquid and of the crystal phases: A G(T,p) = G\,{T,p) — G s (T,p). 
Since it vanishes at the melting point Tm, one can expand AG up to the first order 
of the undercooling AT = Tm - T as 


(3.1) 
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Part I. Ideal Growth Laws 



Figure 3.1: Schematics of the potential surface in configuration space. Crystal phase 
corresponds a stable phase, liquid a metastable phase, and in between is the diffusion 
activation energy. 


In the last equality we used the relation (2.3). 

Under this chemical driving, the crystal grows. For a liquid molecule to be in¬ 
corporated into the crystalline order, it has to change the configuration. But around 
the liquid molecule there is a high density of other molecules and they hinder the 
free motion of the molecule. It can mainly vibrate around its average position with a 
frequency u, which is of the order of that of the lattice vibration. In order for a liquid 
molecule to change the configuration drastically, it has to overcome the energy barrier 
E a of the molecular diffusion, as is shown in Fig.3.1. At a temperature T a molecule 
acquires the energy fluctuation E t \ with a probability proportional to the Boltzmann 
weight exp(-E d /kf/T), where k B is the Boltzmann constant. Therefore, the crystal¬ 
lization rate per unit time is given as u exp(-E^ / k bT) . There is, however, a counter 
effect, namely the melting of a crystal molecule into the liquid state. Since the Gibbs 
free energy per molecule, called chemical potential ft = G/N, is higher in the liquid 
phase than in the crystal phase, /q, > ft s , the rate of melting is smaller than that of 
crystallization by a factor exp(—Aft/AsT). Here A/t is the difference of the chemical 
potentials between the liquid and the crystal phases: Aft = fti,(T,p) — ft s (T, p). By 
the crystallization of one molecule, the solidification front increases its height by a 
molecular height a. Thus the growth rate is given as 


V = au exp 



1 — exp 



(3.2) 


In terms of the liquid viscosity p or the diffusion constant D one can describe the 
formula as 


V = 


k B T 

Qita 1 ri 


1 — exp 



= K 


1 “ cxp (“0) 


(3.3) 


where the Einstein-Stokcs relation [56] 




= D = 


knT 

6nr/a 


(3.4) 




4. Vapor Growth 


7 


is used. Here K — kf{T/6na 2 ri is called the kinetic coefficient. The linear growth law 
(3.2) or (3.3) is called the Wilson-Frenkel formula for the melt growth [196, 63], and 
for the small undercooling AT it is approximated as 

v ~ *0 * KtAt ' (3 - 5) 

where Kt = Kl/ksTTu with l = LjN being the latent heat per molecule. At small 
undercooling AT, the growth rate is proportional to the undercooling. 

The growth rate thus obtained is an ideal one and valid only when the interface 
is atomically rough and the whole surface is accessible for the crystallization. In the 
actual crystal growth the roughness of the interface controls the growth rate. Also 
the derivation of Eq.( 3.2) lacks the consideration of latent heat and of the associated 
temperature increase near the interface. The effect of heat transport is important 
in determining the growth rate and the morphology of the crystal, but that will be 
discussed later in Part III. 

The Wilson-Frenkel formula says that the growth velocity V should drop drasti¬ 
cally at low temperatures, since the liquid viscosity r] increases exponentially. In a 
molecular dynamics simulation of the crystal growth of the simple molecule system 
[39], on the other hand, the growth rate is not limited by the mobility of the atoms in 
the bulk liquid: The kinetic coefficient K in Eq.(3.5) is found to be proportional to 
the temperature. Explanation is given such that the precursor of the crystalline order 
is already formed in the liquid phase near the interface, and the collective motion of 
the liquid facilitates the crystal growth without the activation barrier [137]. Thus, 
general and microscopic consideration on the kinetic coefficient seems worth to be 
studied. 


4 Vapor Growth 


Crystal can grow not only from a melt but also from a gas phase, as shown in Fig.4.1. 
For instance, the snow is a crystallized water from the vapor (Fig.l.2a), whereas the 
ice is the one grown from the liquid water (Fig.l.2b). In a semiconductor industry, 
vapor deposition of a thin film on the substrate is an important technology to fabricate 
materials with a controlled design and new functions. Molecular beam epitaxy (MBE) 
or atomic layer epitaxy (ALE) are among these modern technologies. 

For a gas phase an ideal gas is a good approximation due to its low density. At 
a temperature T and a pressure p, the velocity distribution of a monatomic ideal gas 
follows the well known Maxwell distribution [121] 


PMd ''= (dh) 


3/2 


exp 


(- 


mv 2 \ 

2k b t) 


dv, 


(4.1) 


where v = (v x ,v y ,v z ) is the velocity of an atom and m is its mass. In a unit time, 


/ oo roo rO 

dv x / dv y / dv z n\v,\P(v) 

-oo J—oo J —oo 


rmfcflT 


(4.2) 
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Part I. Ideal Growth Laws 



Figure 4.1: Vapor atoms depositing on a flat surface. 


number of atoms impinge on a unit crystal surface normal to the ^-direction as shown 
in Fig.4.1, and try to crystallize. Here n = p/k B T is the average number density of 
gas atoms. 

Inversely, there are some atoms desorbing from the crystal surface at a finite 
temperature T. The desorption flux is independent of the deposition flux from a 
gas phase, but is a function of the temperature. If the crystal is in the gas with a 
saturation pressure, p eq (T), the deposition rate balances with the desorption rate: 
The desorption flux from the crystal / des is equal to the deposition flux / eq of a gas at 
a saturation pressure p eq {T). Assuming that an atom is cubic with a linear dimension 
a, the net atomic flux in an atomic area a 2 contributes to the crystal growth, and the 
atomic height increases by a. The growth rate is thus written as 


V = a\f - /«,) = 


Q{P Peg) 

\J'hvmk B T' 


(4.3) 


Here ft is the specific volume of a single molecule fl = a 3 . This linear growth law is 
called the Hertz-Knudsen’s formula [77, 116]. 

Since the chemical potential of an ideal gas is written as p G (T,p) = po(T,p 0 ) + 
k B T\n(p/p 0 ) [121], the chemical potential difference of the gas and the crystal is 
written asA/i = p G (T,p) - p s (T) = p G (T,p) - Pa(T,p eq (T)) = ksTl^p/p^), or 
P = PeqC Al ‘ /kH ' - Thus the Hertz-Knudsen formula is represented as 


V = H/ eq 





(4.4) 


The growth rate V is proportional to the driving force Ap for small A p. The kinetic 
coefficient K is now obtained as K = fl/ eq . The relation (4.4) is valid only in the 
ideal situation of a rough surface with the fast transport of materials in the gas phase. 
For a smooth and flat surface, the growth law should be modified as will be described 
in Part III. 
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5 Solution Growth 

It is technically hard to crystallize a material with a high melting temperature, but 
by solving it into some solution crystallization becomes easier at low temperatures. 
For example, NaCl melts at 800°C and melt growth has to be performed at a high 
temperature. If it is solved into hot water, the solution becomes supersaturated by 
evaporation and a large crystal can be grown even at a room temperature. 

In the case of solution growth, there should be a solute atom in front of the inter¬ 
face to crystallize. The probability to find a solute atom at a certain crystallization 
point with a unit volume a 3 is a 3 c for a solution with a concentration c. This atom 
is oscillating around the average position with a frequency t/, and tries crystalliza¬ 
tion. In the solution, however, the solvent molecule makes some chemical bonding 
with the solute molecules. Thus, for the solute molecules to be incorporated into 
the crystalline structure, the solvent molecules have to be desolved from the solute. 
For this desolvation process there is an energy barrier E deso i- Among v trials of so¬ 
lidification, the rate which overcomes the desolvation energy barrier is given by the 
Boltzmann weight exp(—E dasal /k B T). Then the velocity of crystallization is given by 
V cry (c) = au(ca 3 ) exp(-E ieB0 i/k B T). There is an inverse process of melting of crystal 
molecules. Its rate is determined by the temperature, and should balances the crys¬ 
tallization from the solution with the equilibrium concentration c eq : Iwi = Kry^e,,)- 
The net rate of crystallization is 

V = V CTy (c) - V m „i = ua 4 exp ( c - c eq ). (5.1) 

The chemical potential of a dilute solution with a concentration c is expressed as 
= A*>oi(T,Co) + k B T\n(c/co). The chemical potential of the crystal is equal 
to that of a solution with an equilibrium concentration c,, q as tts(T) = /i so i(T, c eq ). 
Then the excess chemical potential of the supersaturated solution is written as 

A/i = /i S oi(T,c) - W (T) = k B T In « k B T (± - 1 ). (5.2) 

The growth rate is expressed in terms of the chemical potential difference A ji as 

vsK fa(i$)- i \ aK & < 5 ' 3 > 

The kinetic coefficient is now defined as K — ua 4 c eq exp(—E deso )/k B T). 

The linear growth law (5.2) is valid only in an ideal case when the interface is rough 
and the whole interface contributes to the growth. Actually by the solidification the 
number of solute molecules in front of the interface decreases, and to compensate this 
solute deficiency, molecules have to be transported from the far end of the solution. 
By including this diffusional material transport, we reach to the different growth law. 
But the detail will be discussed later in Part VI. 
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Exercise: If the transport of materials in the ambient phase is important, it in¬ 
troduces a length scale / which characterizes the effective range for the material to 
be incorporate from the environment. Find in this case the growth formula of an 
[100]-face of a crystal with a unit cell ci\ x <12 x “ 3 - 

Answer: The number of atoms to be incorporated on a unit of [I00]-face with an 
area a 2 a 3 is given by ca 2 a 3 l with a solute concentration c. Since the interface advances 
a height by the crystallization of one atom, the growth rate is described as 

Vi = lv£ler Ei ‘"‘'l kuT (c — c„ q ), (5.4) 

where ft = aua 2 a 3 is the volume of the unit cell [23], 


6 Equilibrium Shape 

So far we considered the growth of a crystal with a flat and infinitely large interface. 
If the interface deforms, the interface area changes and the associated energy change 
has to be considered. Also in the early stage of crystal growth, it starts from a 
small embryo. During growth embryo size increases, and the energy cost at the 
interface changes. Thus the interface energy has effects both on the crystal shape 
and the growth dynamics. In this section we focus on its effect on the crystal shape 
in equilibrium [200, 121, 6], Even though the realization of equilibrium shape is 
difficult due to its long relaxation time, there arc some experiments with very small 
crystals [78, 79, 136], or with a crystal with fast material transport [147]. 

Due to the energy cost at the interface, a crystal nucleus with a finite size and 
shape cannot coexist with a mother phase at a bulk transition point, where the 
chemical potentials of two phases are equal: A n — 0. Finite A/r is necessary to 
compensate the energy increase at the interface and to keep the small nucleus in the 
stationary state. 

Let us assume that z = ( (x,y) describes the phase boundary between the crystal 



Figure 6.1: Crystal surface. 
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and liquid phases (Fig.6.1). Then the free energy gain by the solidification of N s 
atoms is given as 

V An f 

G b = N s (p s -p L ) = --Ap=-— j dxdy((x,y), (6.1) 

where fi and V are the atomic volume and total volume of the crystal cluster, re¬ 
spectively, and the volume integration is taken over the crystal phase (z < ((x, y)). 
The bulk free energy G b is negative for supercooled melt (p s < pi.)- There is an en¬ 
ergy increase due to the formation of an interface. The interface free ene rgy per area 
liPxiPy) is a function of orientation of the interface n = (— p x , -p y , 1)/ Jl + pi +p y 
where p x = d(/dx and p y = dQjdy are the interface gradients. The total surface free 
energy is written as 

G a = J dA-y(p x ,p y ) = j dxdyf(p x ,p v ). (6.2) 

Since the projection of the surface area dA in ^-direction is equal to the area dxdy 
in 22 /-plane as n z dA = dxdy, the free energy density / per unit area in zy-plane in 
Eq.(6.2) is given as 

f(Px,Py ) = 7(Px,Py)iJi+pl +p 2 y (6.3) 

The equilibrium state is determined by the stationarity condition of the total free 
energy G - Gb + G„. The change of the bulk free energy 6\, by the small variation 
of the interface height 8(,(x,y) is written as 


6G b = J <5{(z, y)dxdy s -2A J 6C,dxdy, (6.4) 

where we introduce a parameter A = A/t/2ft. On the other hand, the surface free 
energy varies as 


6G a 



d f d6 ( 

dp x dx + dp y dy 



ffo + f-ep, 

dp x dp y 


J 


— / dxdybC, 


d df d_d£ 
dx dp x + dydp y 


(6.5) 


where partial integration is performed to derive the last equality. By imposing the 
stationary condition 

6G b + 6G S ~ 0 (6.6) 

for arbitrary variation 6((x, y), one gets the Euler-Lagrange equation 


d df d df „ 

2A + — — -f — — 0. 

ox op x dy dpy 


The solution of this equation is 


(6.7) 



and 


21 

dPy 


= -Xy- 


( 6 . 8 ) 
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It shows that variables x and p z are conjugate as well as y and p y are, and f(p x ,p y ) 
behaves as a potential function. The Legendre transformation introduces the new 
potential 

g-f- LT-x - T^-y = / + A xp x + \yp y , (6.9) 

uPx &Py 

which is a function of x and y variables as 


dg = df + d{\xp x ) + d(Xypy) 


A( l^ + | d2/1=AdC 


( 6 . 10 ) 


On integration one gets g(—Xx, —Xy) = AC- The potential g is thus equal to the shape 
C as 

Q f 

C = J = J+XPz + VPy (6.11) 

By rearrangement, it reduces to 


C -xp x - yp y _ 'r(p x ,Py) 

Jl+Pl + p\ a 


( 6 . 12 ) 


Since the normal vector to the interface is given by n = [1 + pi + pl]~ l/2 (-p x , ~p v , 1), 
the point r = (x, y,Q on the interface satisfies the relation 


(r ■ n) — 


7 _ 2fl7(n) 

A Ap 


(6.13) 


The analogous formula holds for a one-dimensional interface of a two-dimensional 



Figure 6.2: Wulff plot and equilibrium shape. 
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Figure 6.3: Local coordinate system; z axis in normal direction, and x\ and x% in the 
local principal directions. 


system. Actual realization of a one-dimensional interface is a step on a crystal surface 
with the step free energy /?(n). The equilibrium shape is determined by the relation 


(r • n) = 


fhP 

Ah' 


(6.14) 


where 1^2 = « 2 is an atomic area. 

Since the l.h.s. of Eq.(6.13), h = (r-n), is the vertical length from the crystal center 
to the interface with orientation n, Eq.(6.13) shows that the length k is proportional 
to the surface tension 7 ( 11 ). By using this relation one can draw the equilibrium shape 
of a crystal as shown in Fig.6.2 (Wulff plot) [200]. Draw a line from a center O to the 
direction n with a radial length proportional to the surface tension 7 ( 11 ), and then 
draw a plane perpendicular to it through the end point: The crystal surface can be 
a part of this plane. By varying the orientation n, the envelope of the perpendicular 
planes gives the equilibrium shape. 

Another way to represent the equilibrium shape is its relation with the interface 
stiffness. In three dimensions, we choose the local coordinate (x\,X 2 ,z) such that z 
axis is in the normal direction of the interface, and xi, x -2 axes in the local principal 
directions with principal radii of curvatures, Ri and R 2 (Fig.6.3). Then the profile of 
a crystal near this point is approximately given as 


T 2 t 2 
_ x 2 

2R\ 2Z?2 


(6.15) 


with Xi = R\di, x-i = RzO 2 . Here the angles 61 and 82 are the deviation of the normal 
direction as shown in Fig.6.3. Then the slopes are p; = dzjdxi = —x,/f?, = — 8 , with 
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i = 1 and 2. Equation (6.7) can be expressed as 


Aft X d df _ Ap y. y. dpj d 2 f 
Q “( dx, Bp, fl Bx, dpidpj 


(6.16) 


Here / = 7 ( 6 * 1 ,0 2 )yl + 6 \ + 0f, and the surface stiffness at x\ = xi = 9\ — 6*2 = 0 is 


d 2 f _ d 2 f g _gl_ 

dp.dpj do^ ' 


(6.17) 


The curvature is 

dpj 

dx, 

Thus Eq.(6.16) reduces to the relation 



Ap _ 711 722 

O ~ /?, « 2 ' 


(6.18) 


(6.19) 


If the surface tension is isotropic, then %j — 7 < 5 p and the crystal nucleus is spherical 


with the critical radius 



( 6 . 20 ) 


obtained from Eq.(6.19) with = R 2 = R c . This radius corresponds to the crit¬ 
ical radius of a three-dimensional nucleus in equilibrium with the undercooling or 
supersaturation, Ap. 


Exercise 1: In two dimensions the interface is defined by the arc length s and the 
angle 6 of the normal vector n from y axis. (Fig. 6 .4) 

(1) From Eq.(6.14), show that the point ( x,y ) on the interface satisfies the relation 

Az(s) = /?( 6 >) sin 0+ /?'(<?) cos 0 (6.21) 

A y{s) = /?( 0 )cos 0 -/?'( 0 )sin 0 ( 6 . 22 ) 

with fd'(9) = dfj/dO and A = Ap/H 2 . 



Figure 6.4: Geometry of a step at a site r with a normal vector n. s is the arc length 
and 9 is the angle of the normal vector. 
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(2) Find the relation 


with the step stiffness 


p- fi+ i£ 


and the radius of curvature p = ds/dd. 


Answer: 

(1) Since n = (sin 0, cos 0) from Fig.6.4, Eq.(6.14) can be written as 

(r • n) = x sin 0 + y cos 0 = ^. 

By differentiating both sides by s, one obtains 

dx dy . 0 ,d8 /?' d8 

— sin0 + — cos0 + (.Tcos0 - y sin0) — = —-p. 
ds ds ds \ ds 

From the geometry shown in Fig.6.4, 


and one gets 


dx = ds cos 0 , dy = -ds sin 0 , 


xcosO — j/sin 0 = 

A 


Solving (6.25) and (6.28), the desired relations (6.21) and (6.22) are obtained. 
(2) Differentiating Eq.(6.28) with s, one gets 

^■cos0 - ^~-sin0 - (a;sin0 + j/cos0)^ = (6.29) 

ds ds ds A ds 

Inserting the values of dx/ds and dy/ds from Eq.(6.27) and by using the relation 
(6.25), one gets the relation 

=cos 2 0 + sin 2 0 = 1. (6.30) 

ds A 

By rearrangement, the desired result (6.23) is obtained. 

Exercise 2: Draw the equilibrium shape of a two-dimensional crystal with a tension 

0(0) = A,(l + ecos40). (6.31) 
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Figure 6.5: Polar plot of the step free energy 0 = 0 O ( 1 + ecos 48) in a dashed curve, 
and the equilibrium shape in continuous curve, (a) e — 0.06, and (b) e — —0.06. 


Answer: From Eqs.(6.21-6.22), the shape is determined by the parametric represen¬ 
tation; 

~x = (1 + ecos40)sin0 - 4esin40cos0 
00 

y = (1 + e cos 40) cos 0 + 4e sin 40 sin#. (6.32) 

Po 

The shape is plotted in Fig.6.5. For e > 0, the direction 8 — 0 corresponds to the 
maximum of the surface free energy 0 and to the minimum of the stiffness 0. The 
equilibrium shape has a pointed corner there because the curvature is large. 


Exercise 3: Often a lattice model is used for a simple description of a crystal shape. 
The space is divided into a regular lattice and each lattice site can be empty or 
be occupied by a crystal atom. For an occupied site an Ising spin variable Si = 1 is 
associated, and for an empty site S< = — 1. There is cohesion between the neighboring 
crystal atoms. If the bond from a crystal atom is not connected to the other atom, 
this broken bond causes an energy cost J. This energetics can be represented by the 
Ising ferromagnetic Hamiltonian [121]: 

W = + (6-33) 


Here z is the coordination number, and z = 4 for a square lattice model. There the 
critical point is exactly known to be at [149] 


fceTc __ 1 

J ~ ln(-/2+ 1) 


1.1346- * • 


(6.34) 
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Figure 6.6: Equilibrium shape of 2D Ising model at various temperatures below T c . 
Normalized temperatures are T/T c — 0.1 (outermost) to 0.9 (innermost) with temper¬ 
ature differences 0.1. 


The exact form of the interface free energy 8 (6) is also known at a temperature T as 
a function of the orientation angle 0 from the (01) direction as [159] 


= | cos0|sinh ‘(cd cos#|) + | sin0|sinh '(a|sin0|) 

1C r> I 


(6.35) 


with 



1 — 6 2 


,1 + '/sin 2 2# + h 2 cos 2 20 


\/i 


(6.36) 


and 


2sinh (J/k B T) 
cosh i (J/k B T)' 


(6.37) 


Draw the equilibrium shape of a square Ising crystal at various temperatures by using 
this interface energy [159, 7[. 

Answer: Using the parametric representation, (6.21) and (6.22), the profile is ob¬ 
tained straightforwardly as shown in Fig.6.6. On increasing the temperature the 
shape becomes isotropic and the size becomes small. 


7 Growth Shape 

According to the general formalism of nonequilibrium thermodynamics, velocity of the 
time evolution of a dissipative system is proportional to the thermodynamic driving 
force given by the gradient of the free energy [142]. Thus the time evolution of the 
interface z — (( x,y,t ) is written by 

d((x,y,t) _ KQ 6G 
dt k B T 6£(x,y,t)' 


(7.1) 
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Figure 7.1: Growth shape by Wulff plot using kinetic coefficient K. 

In equilibrium Eq.(7.1) reduces to the stationarity condition (6.6) for the equilibrium 
shape. By choosing the local coordinate such that the z axis is oriented in the interface 
normal n, one obtains the evolution from Eqs.(6.4-6.5) and (6.16-6.18) as 

+ S>,' (72) 

where we show explicitly the possibility that the kinetic coefficient A' depends on 
the orientation n. The growth law (7.2) is similar to the ideal growth formulae, 
(3.5),(4.4) and (5.3), but the chemical driving force A/r is modified by the curvature 
effect. At the pointed part with positive curvatures, Ri > 0 and A 2 > 0, the growth 
rate decreases by the surface tension effect. This is called the Gibbs-Thomson effect. 


For a sufficiently large crystal the radius R\ and /i 2 arc very large, and the last 
term related to the interface stiffness may be neglected. Without it and assuming 
that Afi is constant all over the interface, one can easily integrate Eq.(7.2) up to the 
time t as 



A'(n)Ap 

k n T 


t, 


(7.3) 


where r is a point on the interface oriented in n direction. Eq.(7.3) is similar to the 
equilibrium shape Eq.(6.13), with the interface tension q(n) being replaced by the 
kinetic coefficient K(n) [46]. Therefore, we can draw the growth shape of a crystal by 
the Wulff plot of A'(n) instead of the 7 (n), as shown in Fig.7.1. Also one can derive the 
relation between the local curvature and the kinetic stiffness A' tJ - = KS tJ + d 2 Kf 86,. 
In an isotropic case, the crystal grows spherically with its radius given by 


r = k B l 


( 7 . 4 ) 
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If the orientation n— dependence of the kinetic coefficient K is so strong such that a 
corner appears in the kinetic Wulff plot, the Gibbs-Thomson term with surface tension 
is no more negligible. The growth shape with both surface tension and kinetics are 
generally treated by Miiller-Krumbhaar et al. [142] and by Uwaha [184]. 



Part II 

Statistical Mechanics of Surface 


So far the ideal growth formulae of the crystal arc discussed in a macroscopic, ther¬ 
modynamic sense. In this and the next parts, we consider the microscopic aspects of 
crystal growth. Since the growth takes place at the crystal surface, it is quite natu¬ 
ral to imagine that the surface structure influences the growth. From a microscopic 
model of the surface configuration, a surface phase transition, a roughening transition 
[43, 194, 195], is found: At low temperatures surface is flat and smooth, whereas at 
high temperatures it is rough in the atomic scale. The roughening phase transition 
induces the faceting transition in the equilibrium crystal shape. 


8 Position of Crystallization: Kink Site 

Before considering the crystallization in an atomic level, we have to identify when or 
in what situations an atom is said to be crystallized. In the gas phase an atom is free 
and makes no bond with other atoms. In a crystal it makes z bonds with neighboring 
atoms to lower the energy. An energy gain for each connected bond is set -2 J. For 
N crystallized atoms, there are totally zN/2 bonds, and the total cohesive energy 
is -JzN, or —Jz per each crystal atom. After melting all the atoms become free 
without any bond connections, and the energy is zero. By assuming that the work 
associated to the volume change is negligible, the enthalpy variation by the melting 
or the latent heat per atom is given by l = zJ. 

If an atom freely moving in the gas phase impinges on the crystal substrate and 
makes a bond of z /2 nearest neighbors with the crystal substrate, there is an energy 
gain of zJ, and the atom can be regarded to be crystallized. The problem then is when 
an atom will acquire z /2 bonds. On a completely flat crystal surface, the adsorbed 
atom cannot make so many bonds. But on a crystal surface, there are various defects 
as shown in Fig.8.1. A flat portion is called a terrace. When the heights of two 
consecutive terraces differ, there is a step between them. A step can be straight as 
it runs in a closed packed direction. It can also bent and change the orientation at 
a kink position. When an atom impinges on a flat terrace, the atom makes only a 
few bonds with the underlying atoms. The bond connection is so weak that the atom 
can meander on the surface. This atom is called the adsorbed atom, or adatom for 
short, and performs surface diffusion. During meandering, the adatom may come in 
contact with an uprising step. Then step provides some additional bonds parallel to 
the terrace, but the number of bonds is still in short of z/ 2. The adatom can slide 
along the step edge by edge diffusion, and may reach to the kink position where there 
are z/ 2 nearest neighbors. At this kink site, the atom is really said to be crystallized. 
The striking feature of the kink site is that it never disappears by crystallization or 
melting; it only slides along the step. For the fast crystallization, therefore, it is 
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Figure 8.1: Surface configuration with terraces bounded by a step with kinks. 

desirable to have many kink sites or steps on the surface. But the surface with many 
steps and kinks is rough with many broken bonds, and it has high surface energy. 
The surface cannot be rough at low temperatures where only the minimum energy 
configuration is allowed. The rough surface is possible only at high temperatures. The 
roughness of the surface is expected to change drastically at a certain temperature. 
Since the surface roughness controls the growth mode of the crystal, the roughening 
phase transition of the surface will be discussed in the next section. 


Exercise: In a lattice gas model for crystal growth, a density variable n* at a 
lattice site i takes a value 1 when the site is occupied by an atom, or it is 0 when the 
site is empty. If a broken chemical bond from a crystal atom costs an energy J, the 
interaction energy is expressed as 


H = J [ui(l - Uj) + (1 - n<)nj] - A/r «<> (8-1) 

<ij> i 

where < ij > means the summation over the nearest neighbor sites pairs, and A/j is 
the chemical potential gain by crystallization. Show that the Hamiltonian (8.1) can 
be transformed to the two-dimensional Ising Hamiltonian 

( 8 . 2 ) 

Z <ij> i 

Here an Ising spin variable = 2n, — 1 takes values +1 or —I, and the field is 

H = A/x/2 [128]. 





(a) (b) 

Figure 9.1: Surface configuration with steps observed by scanning tunneling micro¬ 
scope (STM), (a) Si(lll) [192] and (b) Si(001) faces [178], (Courtesy of M.G. Lagally). 

9 Surface Roughening 

We study the equilibrium configuration of the crystal surface in contact with the 
ambient gas phase, which is realized on the vapor-crystal coexistence curve in equi¬ 
librium phase diagram (Fig.2.2). Along this curve the chemical potentials are the 
same for crystal and gas phases. The temperature variation affects only the surface 
free energy associated to the surface configuration or its roughness. At low tempera¬ 
tures the crystal is polyhedral, as shown in Fig. 1.1, enclosed by atomically flat faces 
with low Miller indices, called facets, since these faces have low energies. But the 
flat face has only a single possible configuration, and has no entropy. If atoms jump 
out from the flat face and sticks on it, there is an entropy gain S connected to the 
possible numbers of positions for atoms to be taken out and to be put on. But there 
is an energy cost E, and the equilibrium configuration at a finite temperature T is 
the one which minimizes the free energy F — E - TS. Thus the surface structure is 
determined by compromising two contributions, E and S [194, 195, 146]. 

9.1 Monte Carlo Simulation 

In order to understand the atomic structure of the crystal surface, it is good to see the 
surface configuration: Seeing is believing. Recent development of various microscopes 
as scanning tunneling microscope (STM) and others make it possible to see the surface 
structure in an atomic level, as shown in Fig.9.1. Theoretically, it is also possible to 
produce atomic configurations of the surface by using various computer simulation 
methods. One method is the Monte Carlo simulation [25, 76], which produces various 
microscopic configurations stochastically under the influence of thermal fluctuation. 
Equilibrium configurations of a surface arc obtained by this method as shown in 
Fig.9.2 [194], 

By way of solid-on-solid (SOS) model for the crystal growth (Appendix A9.1), 




9. Surface Roughening 


23 



Figure 9.2: Surface configurations by Monte Carlo simulation; Parameters are the 
temperatures kaT/2J. [194] 

Monte Carlo simulation procedure is now explained [69], Crystal surface configuration 
is assumed to be described by the integer height variables {h} = {/i(l), • • •, h(N)} 
on a square lattice with N = L 2 lattice sites. With the thermal fluctuation, the 
configuration {h} is realized stochastically with a probability P ({/(}, t) at time t. 

The height h(i) at site i changes to another height h'(i) with a transition prob¬ 
ability w(h(i) —* h'(i)) per time. Then the time variation of the probability P in a 
small time increment At is given by 

P (h( 1) • • • h(i) ■ ■ ■ h(N), t + At) = P (A(l) • • • h(i) ■ ■ ■ h{N), t) 

N 

-J> (h(i) -> h'(i)) AtP (h( 1) • • • h(i) ■ ■ ■ h(N), t) 

i=i 

N 

+ E w (/*'(*) Ki)) AtP {h( 1) • • • h'(i) ■ ■ ■ h(N ), <). (9.1) 

i=l 

The second term represents the probability decrease due to the probability escape 
from {h} to other configurations, and the third term represents the increase due to 
the probability income from other configurations to {h}. In the continuous time limit, 
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A t —> 0, one obtains the master equation of the probability; 


dP((h\ t) N N 

■ ■ V" = -E w (M*) - h'(i))P({h},t) + Y,w(h'(i) -* A(*))• 

6,5 .'=1 .=1 

(9.2) 

Here {/i} f is a configuration different from {h} only by a height at the site i. In order 
to realize thermal equilibrium asymptotically (t —* oo), the Boltzmann distribution 


= Z 1 exp 


mm 

k B T 


(9.3) 


should be a stationary solution of the master equation. Here H({h}) is the Hamilto¬ 
nian defined as 


H({M) = ^E m-hU)\ + &it'£h{i) = E{{h})+ ^1:110), (9.4) 

<i}> i i 

with J being the energy cost of a broken bond, Ap = p c - p s the chemical potential 
difference between the gas and the crystal phases. If the transition probability w 
satisfies the detailed balance condition; 


W (h(i) -> h'(i)) p^m.) = pxn nwd-mh}) ' 

w(h'(i) - h(i)) P^Uh}) CP [ k B T 


(9.5) 


the equilibrium is confirmed. The crystallization rate at the site i is chosen to depend 
only on the chemical potential difference as 


w{h{i) -> h(i) + 1) = exp (^) • (9.6) 

The evaporation rate at the site i is determined from the detailed balance condition 
(9.5) as 


»(/,«) - m - 1) = - 1 - '«(')) = «I> (-^|) . (9.7) 

Here AE = E(h(i) — 1) — E(h(i)) is the variation of the bond energy associated 
to the lowering of the height h{i), and is given by A E = —2 J(n — zj 2) with n 
being the number of nearest neighbor sites with heights lower than h(i) and 2 is 
the coordination number. Surface configurations with corresponding n are shown in 
Fig.9.3 for a one-dimensional example (z — 2). Since the maximum of n is z, the 
maximum evaporation rate is proportional to c‘ J l k " T . Thus by choosing 

^ ~ e A,,/fcflT _|_ e zJ/k„T’ ( 9 - 8 ) 

the transition rate wAt in Eq.(9.1) never exceeds unity. 
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i i i 

n=0 1 2 

AE=2J 0 -2J 


Figure 9.3: Configurations of a one-dimensional interface (z = 2) before evaporation 
of the shaded top atom on the i-th column, n means the number of nearest neighbor 
sites with the height lower than hi, and A E is the variation of the bond energy. 

The real configuration change of the simulation is performed as follows: Select 
a lattice site i randomly, and calculate n and various transition probabilities wAt. 
Then, create a random number u which is uniformly distributed between 0 and 1. If 
u lies between 0 and Ate A ^ fcfiT , the height h(i) is increased by 1. If u lies between 
Ate^ kaT and At(e A ^^ kHl 4- e~^ E ^ keT ), the height h(i) is decreased by 1. Otherwise, 
the height does not change. After trying the height variation N times, every site has 
been tried once on average to change its height, and the one Monte Carlo step (MCS) 

is said to have been performed. The time increases by At. See the sample program 

in Appendix A9.2. 

The average height at the m-th MCS is defined as 

Wro)) = (9.9) 

i=l 

where h(i ; m) represents the height at a site i at m-th MCS. In equilibrium where 
A/« = 0, the height should remain constant, and the thermal average of height is 
obtained by 

1 Mp+ M 

( h ) = T7 E ( Km )>• (9.10) 

m m=M 0 

Other quantities are calculated similarly. 

In some cases, the Monte Carlo procedure can be regarded not merely a math¬ 
ematical algorithm to obtain equilibrium quantities, but also as a physical process 
of relaxation dynamics. For example, the SOS simulation so far explained can be 
regarded as to mimic a crystal growth from a homogeneous ambient gas phase. In 
order to incorporate the kinetic coefficient K, the time increases At/K after one 
MCS. The growth rate (5.2) with additional curvature correction is obtained from 
the simulation by 

v (t) = ^ ((H™ + !)) - (*(«*)» > 

where the time is t = mAt/K. 


(9.11) 
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At very low temperatures (k B T < J) with a small supersaturation, the time 
increment At defined in Eq.(9.8) is so small that configuration change occurs very 
seldom. In order to circumvent the difficulty and alter the configuration efficiently, 
the waiting time method is often used. Details are found in references [69, 188, 88]. 

Simulation results of the equilibrium configuration are shown in Fig.9.2 [TO]. At 
low temperatures, the surface is fiat with few bumps and holes, whereas at high 
temperatures the surface is rough. Thus the change of the surface configuration by 
temperature is obvious. This surface structure change is in fact a phase transition, 
which takes place at a critical temperature. For the real understanding of the phase 
transition, one needs analytical theory on it. 

9.2 Mean Field Theory by Jackson; cr-Parameter 

For an analytical study on the roughening transition, we start from the simplest 
model, where the interface is assumed to consist of a single layer between the semi- 
infinite crystal and the vacuum phases [43, 90], The interface layer is divided into a 
two dimensional lattice cells, and each cell site can be occupied by a crystal atom or 
can be empty. Between the neighboring crystal atoms a cohesive energy of attraction 
is assumed. 

If the layer is almost empty, the layer belongs to the vacuum phase and the crystal 
phase terminates sharply before the interface layer: The crystal has a sharp interface. 
If the interface layer is almost fully occupied, the layer belongs to the crystal, and 
the crystal terminates sharply after the interface layer. Only when the half of the 
layer is empty and the other half is occupied by crystal atoms, the interface layer 
belongs neither to the crystal nor to the vacuum phase. In this case the interface can 
be regarded to be rough. Burton, Cabrera and Frank treated the square lattice case 
exactly [43], because the model is equivalent to the exactly soluble Ising ferromagnet 
as (8.1) or (8.2). Jackson later treated the same model in the mean field approxima¬ 
tion [90], and estimated the roughening temperature Tr , and compared it with the 
melting temperature 7 m • Here I summarize briefly the latter mean-field treatment. 

The interface layer consists on N cell sites, and each cell site has z B nearest neigh¬ 
bor sites in the layer. When there are Ns crystal atoms on the interface layer, the 
energy increase is estimated as follows: Among z, neighbors of each crystal atom, 
(1 — N s /N)z s of them are expected to be empty in the mean field approximation. 
Therefore, the total energy cost E is proportional to the total number of broken 
crystal bonds as 

E = J Za ( 1 - N S /N)N S . (9.12) 

The entropy S is obtained as the logarithm of the number W of the microscopically 
different configurations. Since W is the number of ways to choose Ns sites from N 
lattice points to locate crystal atoms, it is easily obtained as W = N\/\Ns\{N —TV’s)!]. 
By using the Stirling’s formula In N\ « N In N - N, the entropy is obtained as 


S = k B In W « k B [N In N - N s In N s - (N - N s )ln{N - N s )\. 


(9.13) 
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'V=Ns/N 


Figure 9.4: Free energy F as a function of the parameter <P = N s /N at various 
temperatures. Below the roughening temperature Tr, F has double minima, whereas 
above Tr it has a single minimum. 


The dependence of the free energy F = E- TS on the parameter 'P = Ns/N is shown 
in Fig.9.4 at various temperatures. By minimizing the free energy F, we obtain the 
equilibrium value of the parameter ty, or the degree of crystallinity of the interface. 
At a temperature below the roughening temperature 


Tr = 


z e J 
2 fcg ’ 


(9.14) 


there are two minimum points in F-'i diagram, showing that the interface is almost 
empty (if < 1) or is almost full ('P « 1). At high temperatures there is a single 
minimum at the intermediate value $ = 0.5; The interface is rough. 

Since the coupling energy J is related to the latent heat per molecule by l = 
zJ with z being the bulk coordination number, the roughening temperature can be 
expressed as 


2r = 


z s l 
2 zks 


(9.15) 


The comparison of the melting temperature 2 m and the roughening one Tr yields the 
following parameter 

Zul 


a = 


zIcbTm 


(= 

V T M )’ 


(9.16) 


which is readily evaluated by using experimental data. If a < 2 the melting point 
Tm is higher than the roughening temperature Tr and the interface at T M should be 
rough. If a > 2, on the other hand, Tm is below the roughening temperature and 
the interface is flat. The parameter a is called the Jackson’s a parameter [90], and is 
used to classify the material for surface roughening or to guide experiments. 
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This model is equivalent to two dimensional Ising ferromagnet defined in (6.31), 
and the exact solution is known [149], The singularities in thermodynamic quantities 
are well studied, and for example, the specific heat shows a logarithmic divergence. 
But the single layer model of interface cannot correctly describe the roughening tran¬ 
sition. In this model the height difference at two separated points is 0 or 1 by 
definition, and the interface is always flat. If the interface is really rough, heights 
at two distant points should be uncorrelated and fluctuate strongly to diverge at 
infinite separation. Burton, Cabrera and Frank extended the model with multiple 
layers and analyzed it in the Bethe approximation [43], They found no singularity in 
the free energy, meaning the absence of phase transition when the height fluctuation 
increases. The Bethe approximation is better than a simple mean field approxima¬ 
tion, but as for critical phenomena the same qualitative behavior is expected since 
both approximations neglect the effect of fluctuation. For the analysis of the true 
roughening transition, one should include the possibility of arbitrary height differ¬ 
ence at two separated points. The thermal fluctuation influences strongly on the 
phase transition, and one needs sophisticated method as variational [163] or renor¬ 
malization group method [114, 49, 146]. The critical behavior of the phase transition 
turns out to be quite different from that of the single-layer Ising model. However, the 
transition temperature is quantitatively not much different from the two-dimensional 
Ising model, and thus the Jackson’s <*-parameter is used as a convenient criterion for 
surface roughening. 

9.3 Mean Field Argument in terms of Steps 

We now look the surface roughening from another point of view, in terms of steps on 
a crystal surface. For a rough surface there are two- dimensional nuclei, where the 
inside terrace is one atomic height higher or lower than the exterior one, as shown in 
Fig.9.5. For a nucleus with the perimeter step length L, the energy cost is E = JL/a 
due to the broken bonds. Here a is the linear size of the atom. 

As for the entropy, one counts the possibility of step configurations with a fixed 
perimeter length L. A step looks quite similar to a polymer with L/a segments. Each 
segment should be oriented in one of the z s directions of nearest neighbor connection. 
The first segment can be put in one of z s possible orientations. From the second 
segment on, it can be in one of z s — 1 orientations, because it cannot fold back on 
the previous segment. Then the entropy is approximately given by S ss kg ln( 2 s - 
\) L ! a = k B (L/a)\n(z B - 1). Here we have neglected such restrictions that the step 
loop should be closed and cannot cross itself. The free energy is then obtained as 
F = [ J — kgT\n(z s — 1)] (L/a) = f3L, with p being a step free energy per length. 
At a low temperature below the roughening temperature Tr = J/fca ln( 2 s — 1), f3 is 
positive and the state which minimizes the free energy F contains no step (L = 0): 
The surface is flat. At a high temperature (T > Tr) L = oo corresponds to the state 
of the minimum free energy and the surface becomes rough with an infinitely long 
step running. In this case the step-step interaction becomes important and the step 
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Figure 9.5: (a) A two-dimensional nucleus bounded by a step with length L. (b) Top 
view of a step loop. 

free energy is shown to remain zero. 


9.4 Short Summaries of Roughening Transition 

In order to understand the roughening phase transition correctly, the mean field 
approximation is insufficient and the fluctuation effect has to be considered carefully. 
Treatments by the variational method [163] and by the renormalization group theory 
[114, 49, 145, 146] are described in Appendices A9.3-4 in detail, and here we give a 
short summary of the renormalization results. There are also some exactly soluble 
models of the surface structure [187, 93] which shows the roughening transition, as is 
described in Appendix A9.5. The result is also summarized here. 

In the SOS model the interfacial configuration is described by the height variable 
h( r) at a two-dimensional site r. Below the roughening temperature 7\ the interface 
is flat and the height fluctuation is small such that the height difference correlation 

G(r) = ((/i(r + r 0 ) - h(r 0 )) 2 ) (9.17) 


at two separate points remains finite. On the other hand, above T R the interface is 
rough and the height fluctuation is so large that G(r) diverges as r -» oo. According to 
the renormalization group calculation [146], G(r) diverges logarithmically above Tr. 
Below Tr the height fluctuate logarithmically at short distances within the correlation 
length £, but for large distances r > (, the fluctuation remains constant: 


G(r) 


In r, 

In*, 


for r < £ 
for r > £ . 


(9.18) 
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The correlation length diverges at the roughening point T R as 

£ ~ ^P ^Ta'-y l' ( 9 - 19 ) 

Within the scale of correlation length £, the surface is rough, as Eq.(9.18) shows. 
This means that steps with a perimeter length up to ( can be excited by thermal 
fluctuation. By denoting the step free energy per length as ft, the total step free 
energy fj( of a step with a length £ should be of the order of the thermal energy k n T. 
Therefore, the step free energy 0 is given by 

/?~r l ~exp[-^L=]. (9.20) 

The singularity of the step free energy is an essential one, whose derivatives of arbi¬ 
trary order vanish at the transition point T R and there is no divergence in derivatives, 
such as in specific heat. Above T R , the correlation length £ is infinite and thus 0 — 0. 

There are some special models of a crystal surface, which are exactly soluble 
[187, 93], They are mapped onto the inverted F-model, a special case of a six-vertex 
model [129], The step free energy is shown to have the same singularity with that 
predicted by the renormalization group method, Eq.(9.20). 

The Monte Carlo simulations also show that the height difference correlation func¬ 
tion G(r) diverges logarithmically above T R , and it saturates below T R [127, 164]. 


10 Step Fluctuation in Equilibrium 

Above the roughening temperature the crystal surface is rough and there are kinks 
everywhere on the surface. In this case the growth formula of Hertz and Knudscn 
may be valid, if the transport is sufficiently fast. Below Tr the surface is flat and 
there is no thermally excited kinks. Then, how can the crystal grow? 

When an atom impinges on a flat surface from the ambient gas phase as shown 
in Fig.8.1, the adatom moves around the surface before it evaporates back into the 
gas. If an adatom meets other adatoms during surface diffusion, they make bonds 
and form a cluster to gain energy. The cluster is stabler than the isolated adatom. 
Further coalescence of adatoms makes the cluster larger to form a nucleus as shown 
in Fig.9.5. Around this two-dimensional cluster or crystal nucleus there is a step and 
kink positions. Thus the two-dimensional nucleus provides the kink sites for the 
crystal growth. This mechanism of growth is called the two-dimensional nucleation 
and growth. 

There is another mechanism of crystal growth. The molecular arrangement of the 
crystal is not always perfect, but sometimes there arc defects. Among them is the 
screw dislocation, a line defect where the molecular arrangement is shifted in two 
consecutive lattice planes. When the screw dislocation ends up at the crystal surface, 
it produces a step running from the dislocation line (Fig.10.1). If there are many 
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(a) (b) (c) 


Figure 10.1: Spiral growth by a screw dislocation. 



Figure 10.2: Spiral on SiC.(Courtesy of I.Sunagawa.) 

kinks on a step, the meandering adatoms meet the step, are incorporated in the kink 
sites, and the step advances. But since one end of the step is pinned by the screw 
dislocation, the step winds up around the dislocation in a spiral form. This spiral 
step never disappears but grows to supply the kink position. This mode of growth is 
called the spiral growth (Fig.10.1 and 10.2). The nucleation and spiral growths are 
two main mechanisms for the growth of the crystal with flat faces. Growth laws of 
the two modes, nucleation and spiral growths, will be described later. 

First we calculate the kink density on the step. Since the kinks are related to the 
change of orientation of the step, it is related to the step fluctuation, too. Nowadays 
step fluctuation can be observed by various atomic-level microscopes as scanning 
tunneling microscope (STM) [178, 192, 86, 41], atomic force microscope (AFM), 
reflection electron microscope (REM) [4, 125, 126] and so on. STM and REM obser¬ 
vations are shown in Fig.9.1 and Fig.10.3, respectively. From these observations one 
can derive the microscopic information of step stiffness and energy parameters. 

Let us consider a simple model of the step running, on average, in x direction, 
as shown in Fig. 10.4. By denoting a lattice constant o, the step position is written 
as y(i) at x = ia. When the step runs in the next position x — (i + l)a, its height 
y(i + 1) can be the same with the previous one, or can be up or down by a single unit 
and form a kink. The total energy of a step running from x = 0 to x = L is written 
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(a) (b) 

Figure 10.3: (a) A step train and (b) an isolated step observed by reflection electron 
microscope (REM) [4]. 


as 


L L l a 

H = J y ~ + J*E 
a 1 = 1 


\y{j) - vii - i)l 

a 


( 10 . 1 ) 


with J x , J y being the energy cost of a broken bond in x and y directions respec¬ 
tively. Since the height difference (y(i) - y(i - 1)) /a between neighboring positions 
takes only three values, 0 or ±1, one can easily calculate the partition function at a 
temperature T as 


Z = 


+i 


E exp ^ 


n=—1 


Jy T J x | It I 


k„T 


i f/" 




Therefore, the step free energy per length is obtained as 
fcflTln Z 1 




-(j,-t„Tl»[l+2exp (—)]). 


(10.3) 



Figure 10.4: Solid-on-solid model of a step running on average in .x-direction. 
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The density of a kink is also obtained as 


_ 2 exp (— Jx/fcaT) 

1 + 2exp(-J x /k B T)' 


(10.4) 


J x thus corresponds to the kink energy. At a finite temperature, there is a finite 
density (|n|) of kinks on a step, and thus the adsorbed atoms reached to the step 
can be easily crystallized at kink sites. If J x = J y = J, the step free energy density 
Po vanishes at a temperature T R = J/(k B In 2). This is said to be the roughening 
transition temperature of this model. 

The correlation function of the step height difference at a separation x defined as 


x/a 


x( a 


G(x) = {(y(x) - y{ 0)) 2 ) = <(£[y(i) - y(i ~ l)]) 2 } = (d> n <) 2 > 


(10.5) 


i=l 


characterizes the roughness of the step. Since the height difference an, at two neigh¬ 
boring sites i and i — 1 has no correlation with those at other positions, the step 
fluctuation is easily calculated as 


x/a 

G(x) = a 2 ^(n 2 ) = ax{n 2 ). 
1 = 1 


( 10 . 6 ) 


The height correlation diverges in proportional to the separation x as x increases, 
and the step is rough at any finite temperatures. In the present model, (n 2 ) = (|n|) 
and the coefficient in Eq.(10.6) is proportional to the kink density (10.4). 

We can relate the step stiffness 0 of Eq,(6.24) to the step fluctuation. In a coarse¬ 
grained continuum model, the total step free energy is written as 


F,t e p = j ds0 (0(s)) = j dx0 (0(s)) 


1+1 f 

ax 


(10.7) 


Here s denotes the arc length along the step, and 9(s) is the angle between the normal 
vector n and y axis as tan 0 = -dy/dx. See Fig.6.4. As the step is rough, the step free 
energy p is expected to have no singularity in the orientation ^-de pendence. Th us, 
for a small step fluctuation, one can approximate 0- & dy/dx, yj 1 + (dy/dx) 2 « 
1 + \{dy/dx) 2 , 0(0) = 0o + ^(cP0/dO 2 )O 2 . Then the integrand in Eq.(10.7) can be 
expanded up to the second order of the derivative dy/dx as 

F, tep ^ Jdx[0 o + \p(^)% (10.8) 


where 0 = 0+0" is the step stiffness in 0 = 0 direction [60, 61]. From this deformation 
free energy F 9tcp , the thermal expectation of the slope fluctuation is easily calculated 
by the equipartition as ((dy/dx) 2 ) = ksT/0. Since the slope at different positions 
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are uncorrelated, the height difference correlation function G{x) is easily calculated 


as 


G(x) = { 




(10.9) 


Comparison of Eq.10.9) with the previous microscopic result (10.6) indicates that the 
stiffness p is related to the height difference fluctuation as 


P = 


k B T _ k B T 1 + 2cxp(-J x /k B T) 


a{n 2 ) 


2a 


exp (~J x /k B T) 


( 10 . 10 ) 


We choose the zero of the y coordinate to the average position of the step as 
(y(x)) = 0. Then the step fluctuation is characterized by the width w defined by 
w 2 = (y(x) 2 ). In order to calculate w, it is convenient to transform the step fluctuation 
in Fourier modes: 

y(x) = '£y(q)e i “ x , (10.11) 

i 

where q = 2itm/L with the integer m — —oo, • • •, oo. Here we assume the periodic 
boundary condition: y( 0) = y(L). The mode amplitude is obtained by 


V (?) = jl y(x) c <qx dx = !/*(-?)■ 

The step free energy F st( , p is written in terms of the Fourier modes as 

iop = PoL + l ^-Y, ?%(?)| 2 ' 


( 10 . 12 ) 


(10.13) 


The thermal average of the fluctuation of each mode is obtained from cquipartition 
as 

i2\ k B T 


<M?)| 2 > 


LPq 2 ’ 


and the step fluctuation width w is calculated as 


w 2 = 


V^/i ( \|2\ 1 kpT /L\ 1 kf}T 

? <l!,(,)l) = -W^f = Tis(rJ = 


(10.14) 


(10.15) 


Reflecting the roughness of the step, its width w diverges for a long step as L 1/12 . 


Exercise 1: Show that a step pinned at both ends, y( 0) = y(L) = 0, fluctuates twice 
larger than a periodic step, and has the step fluctuation 
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The formula can be used to analyze the fluctuation of an isolated step shown in 
Fig. 10.3b to obtain the step stiffness P [4], 

Answer: Step fluctuation is decomposed into sine-modes as 

y(x) ='£y(q) sin qx (10.17) 

9 


with q = nm/L and m is the positive integer asm = 1,2,• • •,oo. The amplitude y(q) 
is obtained by 

2 r L 

y(q) — — / y(x) sin qxdx. (10.18) 

L Jo 

From the ortho-normal condition 


2 

L 


rL 

/ sin qx sin q'xdx 
Jo 


2 r L 

— / cos qx cos q'xdx 
L Jo 




(10.19) 


the step free energy is written as 


<? 


From the 
as 


■fstep — 4 - 

<? 

equipartition, the thermal average of the amplitu< 

/ { \ 2 \ ^kftT 

m > = 


( 10 . 20 ) 

fluctuation is obtained 

( 10 . 21 ) 


The step width w is then calculated as 

w 2 = ^ j {y{x) 2 )dx = ^(ji(?) 2 > 

A : b T /M 2 ^ _1_ = k B T L 
Lj3 V7T/ m 2 6/3 


( 10 . 22 ) 


Exercise 2: If the tilting field H t is applied on the step model (10.1), the Hamiltonian 
is written as 


h = j- + j y I - 1 )\ _ H y( L / a )-y (°) 

v a x a ‘ a 


(10.23) 


and the step is tilted such that the height at the right end (y(Lja)) differs from the 
fixed height at the left end y( 0). 

(1) Calculate the free energy F as a function of H t . 

(2) Calculate the average slope p = (y{L/a) - y(0))/L — (n) as a function of H t , 
and show that for small H t slope p is proportional to H, as 

p=^, (10.24) 

Pa 

with the step stiffness obtained by Eq.(10.10). The stiff step resists tilting 
against the tilting field. 
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(3) Show that the free energy of a tilted step with a slope p defined by the Legendre 
transformation F(p) = F(H t ) + LH t p/a behaves as 

F(p) = F( 0) + £ Pp 2 (10.25) 

for small slope p. 


Answer: 


(1) By defining n, = (y(i) - y(i - 1)) /a, the height difference at both ends is rep¬ 
resented as y(L/a) - y(0) = a T,'=" n,. The partition function is calculated 
as 


Z = 


p. ( Jy + Jx\n\ - H t n 
2^ exp 1 - - 

n=- I 


k„T 

- ai ("S : )l i+tap ( 


1 L/a 


h \ , H t \ L ' a 

kf)T / C ° h k B T. 


The free energy F is then obtained as 

F = —fcflTln Z = ^(j y - fcflTln 1 + 2e'' /l/ * B7 cosh^]) . 


(10.26) 


(10.27) 


(2) The average of the height difference is calculated as 


(y{L/a) - j/(0)) = -a— = L 


2 e -A/fc H T sinh (Ht/fcflT) 

1 + 2c~ J H k u T cosh(Ht/k B T)' 


The slope is then obtained as 

2 e -J./*»T sin h {H t /k B T) _ H t 2c~ J */ kBT _ H t 
P ~ 1 + 2e- J */ k " T msh(H t /k B T) ~ k B Tl + 2 c^'k B T ~ 


(10.28) 


(10.29) 


(3) Since Lp = —adF(H t )/dHt, the free energy F{p) is calculated as 

F(p) = F(H t )-H t ^^^F(0) + H t F'(0)+ l -HfF"(0) 

- H t [F'(Q) + H t F"( 0)] = F( 0) + ^p 2 . (10.30) 

In the last equality wc used the approximation (10.29) for small p and the 
relation F"{ 0) = -(L/a)dp/dH t = -L//?a 2 . The result (10.30) agrees with the 
free energy of a step with a constant slope p = dy/dx calculated from Eq.(10.8). 
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Exercise 3: Consider the similir lattice model for the step but the step height dif¬ 
ference between the neighboring sites, n = ( y(i + 1) — y(i))/a, can take arbitrary 
integer values from —oo to oo. Then calculate the step free energy and the roughen¬ 
ing temperature. This roughening temperature agrees in fact with an exact transition 
temperature of a square Ising ferromagnet. Also calculate the height difference cor¬ 
relation function (n 2 ) and the step stiffness ft. 

Answer: The partition function is calculated as 


= exp (- 


aksTJ 


£ exp - 


\n\J x 
k B T , 




L/a 


(10.31) 


where the parameters a x and a y are defined as 


&x,y 


k B T' 


Then the step free energy per length 0q is obtained as 

A> a In Z ,,a* 

Sr = "I7J = “>" lnco,h T- 

At the roughening point, /J 0 = 0 and thus 

a y ,R = lncoth^. 

It is straightforward to transform (10.33) in a symmetric form as 


(10.32) 


(10.33) 


(10.34) 


sinha^Rsinhctj,^ = 1. (10.35) 

This is the relation to give the transition temperature of the square Ising ferromagnet 
with the nearest neighbor interactions, J x and J y in x and y directions respectively 
[193]. For an equal coupling, J x = J y = J, the transition temperature is at T R = 
J/k B ln(\/2 + 1), in agreement with Eq.(6.34). 

Fluctuation of a step is calculated as 


<n 2 ) = 


2e~"* 

(1 -e“ a ') 2 


1 

2sinh 2 aj,/2’ 


(10.36) 


0 = 


k B T 

a{n 2 ) 


2kgT , , 2 ^ 
—— smh 

a l 


(10.37) 


and the stiffness as 
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Figure 11.1: (a) Vicinal face and (b) the top view of step configuration. 


11 Step-Step Interaction induced by Entropy 


Consider now a surface with a small inclination, which is called a vicinal face. A 
vicinal face consists of terraces of a low indexed singular face and steps separating 
these terraces, as shown in Fig. 11.1. We consider for simplicity that steps arc 
running in x direction with a mutual separation C. The slope of this vicinal face is 
then \dz/dy\ = a/C, where a is the step height. The surface normal makes an angle 
4> with 2 -axis: 

| 0 | = arctan^ « ^ ( 11 . 1 ) 

for a small tilting. In a unit projected area in .ry-plane, there are 1 /C steps running 
and the surface free energy / is written as 


m = 


7(0) 

COS0 


7o + 


m 
c ' 


( 11 . 2 ) 


where 7 ( 0 ) is the surface free energy per unit area, 70 is that of a flat terrace, and @(<p) 
is the step free energy per length for a vicinal face with a tilting 0. At the absolute 
zero temperature, T = 0, the step runs straight and fia is equal to J y of the energy 
of the broken bond. At a finite temperature, a single step fluctuates and the step 
free energy decreases as Eq.(10.3). For a vicinal face, there are many other steps and 
they cannot cross with each other because of the large energy cost for a configuration 
with overhangs. Thus the fluctuation of a step is limited by the neighboring steps, 
and the entropy is reduced, or the step free energy is enhanced. 

We consider the same step model proposed in the previous section. Since the kink 
formation costs an additional energy J x , the probability that a step meanders plus or 
minus y direction is equal to the Boltzmann weight, 

_ exp (~J x /k R T) 

Pk ~ 1 + 2exp {-J x /k n TY 


(11.3) 
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Then meandering of a step in j/-direction can be regarded as a trajectory of a one 
dimensional Brownian motion with a “time” x/a: At each time step Brownian particle 
meanders to plus or minus in y directions with a probability pk and remains there 
with a probability 1 - 2pk- Then after x/a time steps, particle deviates from the 
original position about (A y/a) 2 ~ ( x/a)pK■ If the step fluctuates a separation t, 
then it collides with the other steps on average. The collision takes place at every 
length of Ax ~ apj/(£/a) 2 . This means in a unit length there are (Aa:) -1 collisions, 
and the step looses the entropy. The step free energy increases by about ksT. Thus 
the step free energy j3(<j>) per unit length can be estimated as 

m - 0o + ck B Tp K X « A + *f T s A + (11.4) 

where flo is the step free energy of an isolated step obtained in the previous section, 
and the second term represents the contribution from the step interaction with c 
being some constant. The coefficient is inversely proportional to the stiffness of 
an isolated step as oc (knT) 2 /fj(G). It reflects the fact that the step interaction is 
brought about by the step fluctuation. 


Exerclse:(Gruber-MulIins model [72]) In order to treat the entropy reduction by 
the step confinement more quantitatively, we consider a step confined between two 
walls. For simplicity we consider every length in unit of lattice constant a here. The 
wall is located at y — 0 and y = 2£, so that the step height y(i) at the site with the 
^-coordinate i can take values only from 1 to 2£ - 1. The system has the interaction 
Hamiltonian (10.1). The partition function is transformed in terms of the transfer 
matrix T as 


Z 


, [ J \ 2T-1 M-l 

exp \~T~t) "■ £ T(y(l),y(2))---T(y(L) } y{l)) 

V 7 »(!)=! y(/,)=l 



Here Tr means the trace of the matrix, and T is the symmetric (2t — 1) x (2£ — 1) 
transfer matrix. 

(1) Show that T is explicitly written as 


T{h,ti) = 


( 1 B 0 
B 1 B 
0 B 1 

V o ' ' 


0 

0 

B 


0 \ 
0 
0 


0 B 1 ) 


( 11 . 6 ) 


by using the abbreviation of the Boltzmann factor a s B = e J */ k i>T 
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(2) If one knows the eigenvalues n-k (k = 1 ~ 2i — 1) of f, the trace is easily 
calculated as 

2 £-\ 

Tr T L = £ d- (11.7) 

«-•=l 


For a long step, L —* oo, only the largest eigenvalue fx i gives the dominant 
contribution, and the free energy F is obtained by 


F = — kfjTliiZ — L[J y - fc B Tln^i]. (11-8) 


By denoting the eigenvectors as ij> = (ip(l),p>(2), • • •, ij>(2i - 1)), explicitly write 
down the eigenvalue equation Tip = nip. 

(3) Since the eigenvalue equation corresponds to the wave equation with a fixed 
boundary condition, the eigenvector ip k is written as 

A-(m) = s in^r (fc = 1 — 2€ — 1). (11.9) 

Find the corresponding eigenvalue fi k . 

(4) Find the step free energy per length fJ(tp) of a vicinal face with an angle (11.1), 
and show that for small (p or large fl{cp) can be expanded as 

3(<p) = 8 0 + fhO 2 . ( 11 . 10 ) 


Answer: 

(2) For m = 1 • • • (2<? — 1), eigenvalue equations are written as 

B [il>(m + 1) - 2 ip(m) + ip(m — 1)] + (2J3 + 1 - = 0, (11-11) 

and ip{0) = i>{2t) = 0. 

(3) From the simple trigonometry, one obtains the eigenvalues 

/j* = 1 + 25 cos ^ (11.12) 


for k = 1 ~ 2£ - 1. 

(4) Step free energy is calculated as 
k H T In Z 


m = 


= Jy — kf)T 111 


1 + 


2 exp ^ 


Jx \ 7T 
-1 cos — 


k B T)^2t 


k B T In 


1 + 2cxp(- 


k B T' 


+ 


2k n T cxp(-J x /k R T) 
1 + 2cxp(—J x /h B T) 


= 00 + 02<P 


c 0 4 ) 

(11.13) 
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where /3 0 is given in Eq.(10.3) and 

it 2 k B T 2exp (~J x /k B T) _ tt \k B T) 2 

8 a l + 2exp (~J x /k B T) 8a 2 0(0) ' 


(11.14) 


In the last equation, the lattice parameter a is explicitly written again. If the 
neighboring steps are allowed to fluctuate, the exact calculation shows [3] that 


^ r\k B T) 

6a 2 /?(0) ‘ 


Exercise 2: 


(1) Show that when the height of the vicinal face deviates z(x,y) from the mean 
position, the surface free energy is described by [146] 


H = Jdxdy +ril( 


with two surface stiffness constants 




and 


Til = 


(itkTf 

0£ o a 2 


(11.16) 


(11.17) 


(11.18) 


Here 0 is the stiffness of the steps running on average in x direction with the 
separation £q i n D direction, a is the jump in the terrace height at the step. 


(2) Show that the width of a single step fluctuation is proportional to the logarithm 
of the system size L as 

w 2 = (h 2 > = ^ In L, 
by using the surface free energy (11.16), 


(11.19) 


Answer: 


(1) When steps are straight and aligned regularly as in Fig.ll.2(a), the vicinal face 
z(x,y) = 0 is tilted from the terrace face by an angle <po — arctan(a/£o)- The 
step fluctuation induces height fluctuation, and we consider two independent 
modes of step fluctuation as follows. 

Let every step deforms locally at a position x in y direction by h(x), while 
keeping the step separation £q constant, as shown in Fig.11.2(b). This step 
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(a) (b) (c) 


Figure 11.2: (a) Flat vicinal face, (b) step bending deformation, and (c) step dilata¬ 
tion. 


deformation corresponds to the height increment z(x, y) = h(x)a/C 0 . For a step 
segment of length dx, the energy increase is given by 


¥ 



( 11 . 20 ) 


with the step stiffness /3. In the interval dy in y direction, there are dy/( Q steps 
and the energy increase in an area element dxdy is obtained as 
2 


A E l = i/? 


where 


dz\ 

2 dy 

i. 

( dz\ 2 

(11.21) 

dx ) 


- ipx 

U) dxdv ' 

7i 


a4> o' 


(11.22) 


The second mode of step deformation is the compression or dilatation of the 
step train, as shown in Fig.ll.2(c). When the straight steps shift in y direction 
and one of the terrace width increases from i 0 to (, the tilt angle of the vicinal 
face there in y direction changes by an angle 8<f> = aictm(dz/dy). This portion 
is inclined from the terrace face by an angle <p = arctan (a/i). Three angles 
cj>, 4>o, t>4> are related as <j>o = 4> + and for a small tilt relations are reduced 
to 


a _ a dz 

1 4 dy 


(11.23) 


Now we calculate the energy change caused by the step compression. As is 
explained in the main text of the section, the step collision induces the energy 
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increase. For a step separation £ 0 , the collision takes place along the step for 
every length of L c determined from the relation i\ « k B TL c /(i. The energy 
increase by each collision is of the thermal energy k B T. In an area dxdy, there 
are dy/t o steps and each step makes dx/L c collisions, and thus the total energy 
increase is proportional to 


, T dxdy 

k,T T,u 


{ k B T) 2 

ml 


dxdy. 


(11.24) 


The precise calculation [3, 161] shows that the energy increase is 7 r 2 /6 times 
larger than the above estimation. Therefore, the terrace width variation leads 
the energy increase of 


A£ii = 


tf 2 (k B T)‘ 


(H) 


dxdy 


_ ir 2 (k B T) 2 


1 1 dz 


1 


(tr k B T) 2 

3 dz 3 /c 

)z\ 2 1 

6/3 

(lady + ( 0 a 2 \c 

>y) a 3 \dy) 


ady) 
dxdy. 


Since the average orientation of the whole vicinal face is fixed, 
r L > dz 


fljy QV 

Jo dy^ = *(Ly) — 2(0) = 0 


dxdy 

(11.25) 

(11.26) 


and the linear term in the slope dz/dy is irrelevant in the free energy (11.25). 
By neglecting the third order term of the slope in (11.25), we get 


with 


AE i 


-\ 2 

dxdy 


1 


7|| = 


1. f dz 
2 711 

{itkTf ( 7 rk B T) 2 


pe Q a 2 


fla 3 


<t>o- 


(11.27) 


(11.28) 


By summing up two contributions, (11.21) and (11.27), we get the Hamiltonian 
(11.16). Note the strong anisotropy of the surface stiffness. For a small tilt 
angle <j> 0 7 j_ diverges whereas 7 || vanishes. This drastic difference of two surface 
stiffness constants is observed recently in He [158]. The product 777 y remains 
finite though: 


- - _ 0 _ ( 7 rfc fl T) 2 

7i7|! a <fio 0a 3 a 4 


(11.29) 


(2) From the free energy (11.16), the Fourier-transformed height z(q x , q y ) has the 
correlation 


<I^(9ii9!/)| 2 ) 


k B T 
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, 2 ' 


(11.30) 
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The height-difference correlation function of the vicinal face is then calculated 
as 

<( 2 (r) - ^(0)) 2 ) = 2 f ^-(\z(q x ,q y )\ 2 )(l - cosqr) = -~=\nr = ~\nr. 

J (2-k y n y7x7|| n 2 

(11.31) 

Since the height z(x, 0) is related to the step deformation h{x) as z(x, 0) = 
h(x)a/Ho, the correlation function of the step deformation is shown to diverge 
for large separation along the step as 

{(h(x) - h( 0)) 2 > = 4{ W ,,0) - *(0,0)) 2 ) - 4 ln *■ ( 1L32 ) 

7r z 

One can also calculated the height fluctuation at a point r for a finite system 
of size L 2 as 

= = <1133) 
This result is translated in the step width as 

w 2 = (ft 2 ) = 4 -■ kB ^~ In L = A in L. (11.34) 

The width of a step in a step train diverges weaker than that of an isolated step 
due to the confinement effect of neighboring steps. 


12 Faceting Transition 

We now study the effect of surface roughening on the crystal shape in equilibrium. 
The equilibrium shape is determined by Andreev’s formulation as Eqs.(6.8) and 
(6.11). Therefore, once the free energy / is known, the equilibrium shape can be 
calculated. 

Below the roughening temperature T R , the steps on a vicinal surface with an 
inclination angle 4) in ^-direction (Fig.11.la) has the step free energy (11.4) or 

m = 00 + M 2 . ( 12 . 1 ) 


Then the surface free energy per xy projected unit area, (11.2), is represented as 

M =lo+ WM = lo + W! + « (12 

cos 0 a a a 


Here we have used the symmetry consideration that the surface with positive and 
negative tilt have the same density of steps and should have the same /. Thus / 
has a cusp singularity at 4> — 0 as shown in Fig.12.la. The vicinal face with a slope 
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(a) (b) 

Figure 12.1: (a) Surface free energy / per xy projected area, and (b) the faceted 
equilibrium shape near y = 0 below the roughening temperature Tr. 

p v = dz/dy = tan 4> « 0 is located at a distance y apart from the center, where y is 
given by Eq. ( 6 . 8 ) as 

- ^ (A) + 3 ^ 2 ) sgn(d>). (12.3) 

Here sgn(d>) represents the sign of the slope 4>. Since A, fh and fh are positive for 
T < T r , the slope is positive for y < 0 and negative for y > 0. Let us consider the 
region of y > 0 and <p < 0. (See Fig. 12 . 1 .) Then Eq.(12,3) is solved for <f> as 

(12.4) 

for y larger than yo s Po/\a. The height z is determined from Eq.(6.11) as 

A z = f + \yp y - 70 - 2a' 1 P 2 \<t>\ 3 ■ (12.5) 

The crystal shape is then written as 

for y > yo- Here zo = 7 o/A. For y < -y 0 , <j> is positive and the profile is similar to 
Eq.(12.6) with only a replacement of y to -y. Therefore, between -y 0 and y 0 the 
shape consists of a facet face, and then it continues to the smooth curved faces with 
the exponent 3/2, as shown in Fig.12.lb [160]. The facet size 2y 0 is proportional to 
the step free energy density, ftp. Nea r the roughening point Tr , pp and thus y 0 vanish 
smoothly as y 0 ~ exp [-c/s/Tr - T], This extinction of the facet face is called the 
faceting phase transition, and is equivalent to the roughening transition of the face. 

There are experimental studies on faceting transitions of 4 He crystals in superfiuid 
liquid. The transition temperature is determined to be T R = 1.3K for (0001) face. The 
roughening temperatures of other faces are also determined [ 68 ]. 



Part III 

Kinetics-Limited Growth 


In this part, the growth laws governed by the surface kinetics on a fiat crystal surfaces 
are studied. When the crystal grows with a flat surface, the growth is controlled by 
the incorporation of atoms to steps and kinks provided by the two-dimensional nuclei 
[12] or screw dislocations [62]. Since the step density depends on the growth condition, 
the growth law is different from the ideal ones. 


13 Step Advancement by Surface Diffusion 

We consider the crystal growth from the vapor phase at a relatively low temperature 
where a crystal surface is flat. Atoms deposited on the surface diffuse around on the 
surface before they reach steps. The advancement of a step in the surface diffusion 
field of adsorbed atoms is studied systematically by Burton, Cabrera and Frank [43]. 

We assume that the adsorption site on a crystal surface forms a square lattice with 
a lattice constant a. From the ambient gas phase / atoms deposit on a surface per unit 
area in a unit time interval. Therefore, on a lattice site (x,y) with atomic area Si 2 = 
a 2 , /Sl 2 At atoms are deposited in a time interval A t. The adsorbed atom vibrates 
with a frequency v of the lattice vibration around the averaged adsorption site, and 
tries evaporation or surface diffusion. Due to the coupling with the underlying 
crystal atoms, the adatom has an energy lower than in the free state by £?ad, and the 
evaporation takes place with the probability cxp(-E M \/knT) during v trials per unit 
time. The lifetime until the evaporation is estimated as 

r = j/- 1 cxpiE^/kpT). (13.1) 

By denoting the adatom density as c[x,y), the probability that a lattice site (x, y) is 
occupied by an adatom at a time t is c(x, y; f)fl 2 . The number of adatoms evaporating 
within a time interval A t. is given by c(x, y\ f)f2 2 ■ At/r. 

On the other hand, an adatom hops to one of the neighboring sites during its 
lifetime. For the random walk of the adatom from one adsorption site to the other 
it has to cross over the activation energy or energy barrier E k a, and the probability 
of hopping to the neighboring site is w — ucxp(—E aii /knT). The number of atoms 
coming from the four neighboring sites during the time At is wAfn 2 [c(.T + a, y;t) + 
c(x — a, y\ t) + c(x, y + a\t) + c(x, y — a; <)], while the number of atoms escaping from 
the site (x,y) is wAtQ. 2 Ac(x,y\t). Combining all these processes, the. variation of the 
number of averaged atoms at the site (x,y) during the time interval At is written as 

c(x, y\t + A t)Q 2 = c{x, y\ <)ff 2 + /fl 2 A< - c(x, y ; f )0 2 ^ 

+ w0 2 At [c(x + a,iyt) + c.(x — a,y,t) + c{x, y + a\t) + c(x,y-a;t)-4c{x,y,t)}. (13.2) 
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After dividing Eq.(13.2) by fi 2 At, take the continuum limit At —> 0 and a —> 0 and 
leave up to the second derivative of the space. Then one obtains the diffusion equation 
including the contribution of deposition flux / and the desorption lifetime as 


dc(x,y\t) 

dt 


D„V 2 c+ f 

T 


(13.3) 


Here 

D B = wa 2 = va 2 exp(—E a< i/kBT) (13.4) 

is the surface diffusivity. 

How far can the adsorbed atom meander while it resides on the surface? During 
the time t, an adatom makes wt jumps on average. Since it can jump to four nearest 
neighbors with the same probability, the average position (r (t)) remains the same with 
the original one, (r(0)). But the deviation {(r (t) — r(0)) 2 ) increases in proportional 
to the number of jumps. Since it hops a distance a in a single jump, the deviation 
is written as ((r(t) — r(0)) 2 ) = wta 2 = D s t. Therefore, during its lifetime r on the 
surface, an adatom meanders in the range 


Xs —■ 



(13.5) 


which is called mean displacement of the adsorbed atom or the surface diffusion 
length. 

During the surface diffusion an adatom reaches a step and then crystallizes by 
being incorporated in the kink sites. Let us consider the advancement of an isolated 
step running on average in the ^-direction. Usually the step advances so slowly 
compared to the relaxation of the adatom density c that the stationary distribution 
of c can be realized: The time derivative in Eq.(13.3) is disregarded. If the step is 
straight in the ^-direction, the concentration variation takes place only in y-direction, 
and the diffusion equation is simplified as 


„ d?c _c-c oa 
°dy 2 ~ t 


(13.6) 


where Cx, 
as 


= fr is the concentration far from the step. The solution is easily obtained 


C C 0 o 



(13.7) 


with the surface diffusion length x a . Since the step is rough at finite temperatures as 
shown in section 10, the kinetics of adatom incorporation is expected to be very fast. 
In the extreme case, the step can grow with only an infinitesimal driving as c(0) = c eq . 
From this local equilibrium boundary condition at the step, the integration constant 
A is determined to be A = c eq - Coo, and the diffusion field is obtained as 



c(“T,2/,t) — Coo (Cqo c^) exp 


X. 


(13.8) 
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c 



Figure 13.1: Concentration distribution around the step at y — 0. 
as is depicted in Fig.13.1. 

The net number of atoms impinging to the step at y = 0 from the right (y > 0) 
is wO .2 [c(x % a\ t) — c(a - ,0;t)] = wa 3 dc(x,y)/dy | y=+0 , and that from the left (y < 0) 
is uifl 2 [ c(x , — a;f) - c(:r,0;<)] = -wa 3 dc(x,y)/dy | 0 . Here y = +0(-0) means in 
front (in the back ) of the step. Since these deposited adatoms crystallize at a kink 
site on a step, and step advances a distance a by each crystallization, the advancement 
rate of a straight step is obtained as 



By inserting the density profile (13.8), the step velocity is obtained as 

2D 

vo = (c a o - c oq )n 2 —- = (/ - /«i)n 2 ■ 2£ s . (i3.io) 

Xu 

Here / eq = c eq /r is the equilibrium deposition flux. The result (13.10) can be inter¬ 
preted as follows: The adsorbed atoms landed on a surface within a distance of their 
mean displacement x s from a step will eventually meet the step during their lifetime, 
and be incorporated to it. Therefore, / x 2x s a atoms impinged in a surface region 
of a step length a within a distance x s from the step will crystallize at a step in a 
unit time. As an inverse process, / cq x 2x*a atoms melt back from the step onto the 
surface. Since each crystallized atom pushes the step ahead by a distance a , the 
advance rate is given as v 0 — (J2x K a — f„ t 2x a a) x a in agreement with the result of 
the detailed calculation (13.10). 

So far we considered an isolated step on a singular flat surface. When the surface 
is tilted with a small angle, there are many steps arranged periodically: The vicinal 
surface consists of a step train. When the separation £ between periodic steps is 
smaller than the diffusion length the diffusion field of neighboring steps overlaps 
and competes. A single step gets a contribution only from the small range around 
itself with a width L Therefore the advance velocity will decrease to v = (/ — 
f eq )a 2 £ = Vo x {£/2x s ). Of course, when the step separation £ is wider than the 
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surface diffusion length x„ each step advances with a velocity v 0 as if it is an isolated 
single step. In a more precise treatment, the solution of the diffusion equation (13.3) 
is written as c(y) - c 0 „ = Aexp{-y/x s ) + Bexp(y/x s ) within a region 0 < y < t. By 
taking into account the local equilibrium boundary conditions at neighboring steps, 
c(y = 0) = c(y = i) = Ce q , the integration constants A and B are determined, and 
the concentration is obtained as 


c{x, y; t) = Coo 


_ . cosh \(2y - £)/2x s ] 
l Coo c eq cosh(f/2r s ) 


(13.11) 


From the material conservation (13.9), the advance velocity of a step is given as 


v = (/ - / eq )a 2 2z fi tanh = fotanh ■ (13.12) 

If the step separation £ is much larger than the diffusion length x a , the advance rate v is 
equal to the isolated step v ti in (13.10). If t is smaller than x„ then the approximation 
tanh(£/2x s ) w (/lx, holds and the velocity is approximated as v « (/ - /,.,,)o, 2 as 
already stated. 


Exercise: Derive Eq.(13.12). 

Answer: By imposing a local equilibrium boundary conditions at y = 
one gets 

0 and y — 

c(0) - = A + B = c eq — Coo 

c(() - Coo = Ae~ l/ *‘ + Be t/x ‘ = c cq - Coo, 

(13.13) 

and the integral constants are obtained as 




A (c eq c °°)2cosh(f/2:r s ) 
an d 

e -l/2x. 

(c eq Coo ^2cosh(f/2r B )’ 

(13.14) 

(13.15) 


By putting these results, one gets (13.11). From Eq.(13.9), we get the velocity 
(13.12). 


14 Advancement of a Circular Step 

When a two-dimensional crystal nucleus grows in a diffusion field, the encircling step 
expands and the total step free energy increases. In order for the nucleus to achieve 
further expansion, this energy cost has to be overcome. How will then the step 
advancement rate be modified from the straight step? 
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The chemical potential of the bulk crystal ps agrees with that of the adsorbed 
atoms with the density c«, in equilibrium: p s = Pad(c eq ). If a circular nucleus with a 
radius p is in equilibrium with adatoms, the radius corresponds to the critical one since 
it does not grow or shrink. Then the chemical potential Pad(c eq (p)) for adatom with 
density c K ^{p) should satisfy the critical condition (6.14): A/t = /taxi(c,. q (p)) — /is = 
Here the step free energy fJ is assumed isotropic. Since the ideal solution 
approximation (5.2) is applicable for a dilute adatom system, the chemical potential 
difference is written in terms of adatom density as 

c eq (p) = c eq cxp (jJL) = c eq exp (~^j » c cq (l + ■ (14.1) 

The last approximation is valid for a step with a small curvature k — 1/ p. The curved 
step needs a higher equilibrium concentration than the straight step does in order to 
remain in equilibrium. The effect is called the Gibbs-Thomson effect. 

The adatom density far from the nucleus is determined by the deposition flux 
as Coo = /n. When the density c,*, is higher than the equilibrium value Ceq(/>), the 
circular nucleus with a radius p grows. The density distribution around the circular 
nucleus is expected to be symmetric and to depend only on the radius r. The diffusion 
equation is simplified as 


dc(r,t) /d 2 c 1 dc \ c-Coo 
dt a \dr 2 rdr) r 


(14.2) 


For the slow growth, the stationary approximation dc/dt = 0 is expected to hold. 
The exact solution of c is described in Appendix A14. Here I derive the circular step 
advancement in approximation. 

The radius of the nucleus p is assumed large. Since the step advancement is 
controlled by the diffusion field around the step, the radial variable r in Eq.(14.2) is 
also large, and the second term, r~ x dc/dr , may be neglected compared to the first 
term, d 2 c/dr 2 . The stationary density distribution is obtained as 


c(r) - c x = 


yl +e -(r-e)A. 

J 4_ e ( r -p)A. 


for r > p 
for r < p 


(14.3) 


with x H = \JD h t being the surface diffusion length defined in Eq.(13.5). The adatom 
density relaxes to c x far outside (r » p) and far inside (r <c p) of the cluster. By 
imposing the local equilibrium condition c{p) = c eq (p) at the circular step, the integral 
constant A ± is determined as A± = c eq (p) - Coo- Then the step advancement velocity 
in the radial direction is obtained as 
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where v 0 = 2D 8 Q 2 (c (x - c^/x, is the velocity of the straight step (13.10), and 

P&2 (3^2 (' 14 'i'l 

Pc k B T( Coo /c eq -1) k B T(f/f eq -1) 1 j 

is the critical radius corresponding to the deposition flux / or Eq.(6.14). Therefore, 
the nucleus with critical radius p c does not grow or shrink, the one with radius p > p c 
grows, and the one with p < p c shrinks. The correct calculation in Appendix A14 
leads to the same conclusion. 

15 Growth Rate by Spiral Growth Mechanism 

The growth rate V of a crystal by spiral growth mechanism is now ready to be 
calculated. The step is running out from the center of the screw dislocation. As 
the crystallization proceeds, the initially straight step turns round the dislocation, 
and the curvature at the center increases, as shown in Fig.10.1. But it should be 
smaller than the critical one, k c = l/p c . Otherwise, the crystal melts back due to 
the Gibbs-Thomson effect. Therefore, the curvature at the center is just the critical 
one in the steady state. The simplest form of the spiral is that of Archimedes, 
represented in the two-dimensional polar coordinate (r, (j>) as r = a<p with a to be 
determined. As is explained in the last paragraph of Appendix A15, the curvature 
at 4> -* 0 is calculated to be k = 2/a. Therefore, the parameter a should be 2 p c in 
order to satisfy the boundary condition. The step separation i for large r is given as 
the radius difference by one turn, or A<p = 2v rotation as 

l — a • 2?r = 47rp c . (15.1) 

The more sophisticated calculation is described in Appendix A15 [44]. It shows 
that the asymptotic step separation for large r is 


( « 19p c . (15.2) 

In both cases, step separation t increases in inversely proportional to the driving force, 
( ~ Ap~*. Since the advancement velocity v{i) of a step train is given by Eq.(13.2), 
the time T required for the step to wind up one turn is given by During this 

time T the crystal surface grows a height a in the normal direction, and its growth 
rate of the crystal is then given as 

V = v(i)j - a\f - / eq )^tanh ^ , (15.3) 

and the dependence on the driving force Ap is shown in Fig. 15.1. 

For a low vapor pressure p, the step separation t is larger than the surface diffusion 
length x„ and tanh(^/2x 8 ) « 1. The growth rate is approximately given by V = 
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/// cq - 1 = exp(Ap/k B T) - 1 


Figure 15.1: Growth velocity V versus driving force a = exp(A p/ksT) — 1 for the 
spiral growth. 

a 3 (f - / eq ) 2x s /( proportional to the square of the driving force (p — p M ,) 2 ~ Ap 2 . 
Thus, at low vapor pressure growth follows the parabolic law, different from the 
linear growth law of Hertz-Knudscn. When the vapor pressure p increases, the driving 
force for the crystal growth Ap becomes large, and the critical radius p c and the step 
separation ( becomes small. For t < x s , the approximation tanh(f/2ar 8 ) « f/2x s 
holds, and the growth rate reduces to the Hertz-Knudsen formula given in Eq.(4.4): 
Fhk = a?(f - / cq ) = a 3 (p - p eq )/\/2nmk B T: There are many steps and kinks on the 
surface, and the growth rate V reduces to that of a rough surface. 

16 Two-Dimensional Nucleation and Growth 

Spiral mechanism governs the growth of an imperfect crystal under a small driving 
force. When the supersaturation increases, nucleation becomes important. Also for 
a perfect crystal without dislocation, nucleation controls the crystal growth. In this 
section the growth rate of a crystal by nucleation mechanism is studied [12, 80, 69]. 
The formation of a two-dimensional (2D) nucleus increases the free energy as 

- A//~~ + 2np(i = —Ap ■ n + 0 \J iTrnll 2 — G(n), (16.1) 

where p is the radius of the 2D nucleus, A p = pa — p$ the chemical potential gain 
by the crystallization, the area of a nucleating atom, f) the step free energy per 
length and assumed to be isotropic here. There are n = itp 1 ji\-i atoms in the circular 
nucleus of radius p. G(n) has a maximum as shown in Fig.16.1 at the critical nucleus 
size n c : 

?rpc _ 

n 2 (Ap) 2 ' 


n, 


( 16 . 2 ) 
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Figure 16.1: Free energy cost G(n) of the formation of the two-dimensional circular 
nucleus with the size n. 



Figure 16.2: Variation of the cluster size n by capturing or by evaporating a single 
adatom. 

The maximum value of the free energy barrier is 

G c = G(rtc) = A/m c = (16.3) 

Let us consider the distribution C(n,t) of the cluster of a size n at time t. In 
equilibrium thermal fluctuation allows the formation of nucleus of a cluster size n 
with the Boltzmann weight C°(n) = c,^ exp(—G(n)/fc B T), where c M is the density 
of isolated adatoms: = fr. Note that the equilibrium density of infinitely large 

nucleus is infinite. During the steady growth of a crystal, in fact, a large nucleus 
disappears by the completion of a layer and the density of an infinitely large cluster 
in a steady state should be zero C(n —> oo) = 0. 

The size variation of the nucleus is assumed to take place by capturing or evapo¬ 
rating a single adatom as shown in Fig.16.2. The rate equation is written as 

— w+(n- l)C(n- l) + W-(n + l)C(n + l) — w + (n)C(n) -W-(n)C(n), (16.4) 

where w+(n) is the rate of capture of a single adatom by the nucleus with the size 
n and w-(n) is the rate of evaporation of an atom from the nucleus of size n. The 
detailed balance condition 

w + (n)C°(n) = W-(n + 1 )C°(n + 1) 


(16.5) 
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confirms the approach to equilibrium. Then the ratio of the capture and evaporation 
rates is obtained as 


w+(n) _ C°(n + 1) _ exp[-G(n + 1 )/k B T] 

w-(n +1) C°(n) exp[-G(n)/A:/)T] 


(16.6) 


In a steady state, dC(n)/dt = 0 but there is a net flow J of the cluster growth sup¬ 
ported by the constant supply of particles from the smallest size and their annihilation 
at the largest size. J is then independent of the cluster size as 


J = w + (n)C(n) — w-(n + l)C(n + 1) = w + (n)C°(n) 


C(n) 

C°(n) 


C(n + 1) 
C°(n + 1) 


(16.7) 


For n = 1 there arc plenty of isolated adatoms and one assumes that G(l) = G°(l). 
As for the boundary condition at n — oo, the infinitely large nucleus disappears by 
the completion of a monatomic layer and thus C(oo) = 0. Then by using the relation 


£ 


j 


= £ 


C(n) C(n + 1) 


„=t w+(n)C°(n) LC°(rz) C°(n + 1) 
the nucleation flux J is obtained as 

1 


C(l) C(oo) 
G°(l) G°(oo) 


= 1 , 


J - 


.S w+(n)c°(n) 


-l 


(16.8) 


(16.9) 


Since G°(n) _1 = c” 1 exp[G(n)/fe fl T] has the maximum at the critical nucleus size n c of 
Eq.(16.2) and w+(n) varies slowly in proportional to the surface area as w+(n) ~ n 1 ^ 2 , 
the dominant contribution to the summation in Eq.(16.9) comes near n c . For small 
Ap, n c is so large that the summation can be replaced by the integration as 


£ 


i 

w + (n)C°{n) 


_ 1 _ f dnexp (Gc+jG^in-n^l 

w+(nc)Cco J P V k B T 

1 f G c \ l2nk B T 

w + (n c )c x eXP \k B T)V |G (2 >| ' 


(16.10) 


Here G c is the maximum of the free energy cost given in Eq.(16.3) and G® is the 
second derivative of the activation free energy at n c as 



ln c 


Ap 3 

2tt fl 2 /? 2 ' 


(16.11) 


Therefore, the nucleation rate is written as 


J « Zw + (n c )Coo exp (~G c /k B T) 


(16.12) 
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with the Zeldovich factor 

JW ( A„» Y'* / Ap 

V 2rt,T V ^Hhki_:r;p) Urt.TnJ ' 

Since the adatom capture by the nucleus takes place at the periphery of the nucleus, 
w + (n c ) is proportional to the perimeter length as [80, 69] 


w + (n c ) = 



vo,+ 


2irp c 

Q2 


= 2xJhf 


27T/? 
A p ’ 


(16.14) 


where %+ = 2r g fl2/ is the step advance rate by atom incorporation in Eq.(13.10). If 
there is a direct incorporation of atoms from the ambient gas phase, this contribution 
should be included in wq,+ [69]. Thus 


7 2ar a Q 2 /Ccc / Am\ 1/2 „ „ ( G c \ 

US?) (1615) 

When the nucleation rate J is large, nucleation starts at various places on a crystal 
surface and those nuclei spread with the advance velocity v , coalesce and complete 
a single layer. This mode of layer growth is called the multinucleation growth. We 
now calculate the normal growth rate V of the crystal in this situation. On a flat 
surface it takes a time T for the completion of a single layer growth. During this time 
interval T, nucleation occurs with the rate J per unit area and time. At a time t 
after the nucleation started, its radius is vt and occupies the area n(vt) 2 . Therefore, 
during T the area scanned by nuclei is calculated as 

fT 'Jr 1 

/ (SJ)rr(vt) 2 dt = S-^-v 2 T 3 , (16.16) 

J 0 <3 


where S is the total area. When the scanned area (16.16) coincides with the total 
area 5, the single layer growth is completed, and thus the time T is determined by 

^-v 2 T 3 = 1. (16.17) 

After this time T the crystal surface moves up a distance a, and the normal growth 
rate by the two-dimensional nucleation mechanism is given by 



?r Jv 2 ^j^ 


Since the advance rate of the step is given in Eq.(13.1Q) or 
v = 2x„n 2 f cq [exp (0) - 1 ], 


(16.18) 


(16.19) 
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///cq - 1 = expOty! /k B T) - 1 

Figure 16.3: Schematic diagram of the growth velocity V by the two-dimensional 
nucleation. 

the normal growth rate is 



The last approximation is obtained for small supersaturation A/< where c 0 c « c eq , 
and the result is schematically shown in Fig.16.3. Since V is proportional to the 
exponential of A/i~ l , it is very small for small Ayr. The rate V becomes observable 
only when G c ( Ap) is of the order or smaller than the thermal energy kgT; 

A/t > A/t c = (16.21) 

For large enough A/t, the nucleation barrier vanishes, and the growth rate is V ~ 
(A/t) l,,6 exp(4A/t/3fccT), larger than the linear growth law V <x e A, ^ kBT — 1. This 
result is incorrect. For large A/t, the critical nuclear size n c becomes so small that 
n c < 1, and the crystal surface becomes kinetieally rough. Then the growth rate 
should be that of the rough surface or the Hertz and Knudsen formula (3.4). 

At the thermal roughening temperature Tr and above, the step free energy [i 
vanishes, and the energy barrier G c disappears. Thus the growth mode varies from 
the nucleation controlled exponential type, V ~ exp(-C/A/t) below the roughening 
transition temperature Tr, to the normal growth rate V ~ Ap above Tr. This is 
used for the criterion to determine the roughening temperature in the Monte Carlo 
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Figure 16.4: Growth rate R versus chemical potential difference of the SOS model. 
Numbers aside the curves represent inverse temperatures, 6 J/k B T. k + here is the 
deposition rate, and d is the atomic height. The dashed lines are the growth rates 
calculated from the two-dimensional nucleation model [70] 




(a) (b) 

Figure 16.5: Growth rate of He crystal near and below the roughening temperature 
T r = 1 .23K (a) in the normal and (b) in the semi-logarithmic plots [198]. H here is 
proportional to the chemical potential difference A/i. 


simulation as is shown in Fig.16.4 [70], and in the experiment on He as shown in 
Fig.16.5 [198], 
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Exercise: When the nucleation rate J is very slow or the step advancement velocity 
v is very fast, there is only a single nucleus during the layer growth. Then the normal 
growth rate V differs from multinucleation case. Calculate the normal growth rate V 
of a surface with an area S. This growth mechanism is called the single nucleation 
growth. 

Answer: When a single nucleus is created, it spreads over the whole surface S 
quickly before further nuclei are formed on the same level. The time T necessary for 
the completion of the normal growth of a height a satisfies the relation 


JST = 1. 


(16.22) 


The growth rate is obtained as 


V = — = aSJ = ix^fSc^ 


/AM 1/2 

7T/? 2 fl 2 

te) exp 

Ap/CflT 


(16.23) 


If one plots In V versus 1/A /j, the slope by a single nucleation mechanism is three 
times larger than that by a multiple nucleation mechanism. This crossover is said to 
be observed in polymer crystallization [81,179]. 


17 Asymmetry in Attachment Kinetics 

In the treatment of adatom diffusion, Burton, Cabrera and Frank (BCF) [43] assumed 
the fast kinetics at the step such that the local equilibrium is realized there. Also 
they assumed the symmetry in the attachment kinetics of adatoms from the upper 
and lower terraces at the step. However, the incorporation rate of an adatom from 
an upper terrace needs not to be same with that from a lower terrace as shown in 
Fig.17.1 [172, 190]: From the upper terrace an adatom has to break many chemical 
bonds with the underlying substrate atoms when it crosses over the step down to 
the crystallization position. From the lower terrace, on the other hand, an adatom 
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Figure 17.1: Asymmetry in the incorporation of an adatom (a) from the upper and 
(b) from the lower terraces, called Schwoebel effect, (c) The potential energy profile 
of the adsorbed atom. 
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can simply make additional bonds with the step atoms before it reaches kink sites. 
Additional energy barrier for the crystallization makes the kinetic coefficient from 
the upper terrace A'_ smaller than that from the lower terrace K + . This asymmetry 
in attachment kinetics at the step is first studied by Schwoebel and Shipsey and is 
called the Schwoebel effect [172]. 


17.1 Schwoebel Effect 


We consider the step down configuration, where the step is running on average in x 
direction at y = 0 and the terrace in front at y > 0 is lower than the terrace in the 
back y < 0. Since the step is thermally rough, the step advance rates, v + and 
by the adatom incorporation from the upper and the lower terraces respectively are 
linearly proportional to the supersaturation as 


v± ~ K± 


( PCh \ 

( C± Ceq *■> k B T K ) ' 


(17.1) 


Here the Gibbs-Thomson effect of curvature is included. The Schwoebel effect means 
that the kinetic coefficients are different: K + ^ AT_. 

We now restrict ourselves in the extreme and the simple case that 


AL = 0 and K + = oo. 


(17.2) 


In this case, there is no crystallization from the upper back terrace. Since the crystal¬ 
lization takes place by the atom incorporation only from the lower terrace, it is called 
a one-sided model. Furthermore, the kinetics from the lower front terrace is assumed 
extremely fast such that the local equilibrium is realized: The adatom density in front 
of the step c + is equal to the equilibrium value with the Gibbs-Thomson effect 


— Ceq 



(17.3) 


As in the BCF model, the atoms are deposited on the crystal surface with a flux /, 
the adsorbed atoms then diffuse on the surface with a surface diffusion constant D a , 
and then evaporates back into an ambient vapor after a life time r. For a straight 
step its advance velocity v 0 is calculated to be 


Vo = {f ~ /eqKfl 2 , (17.4) 

where /«, = c^/t is the equilibrium deposition flux, x, = y/D a r is the surface diffusion 
length and fi 2 the atomic area. The velocity vg is half of that given in (13.10) in BCF 
theory, since there is a contribution only from the front terrace. Eq.(17.4) shows that 
the step incorporates atoms deposited in the range x„ in front of the step. 

We now consider the stability of the straight step [8]. If the step is pushed forward 
at some part by fluctuation, the region to incorporate adatoms expands radially as 
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Figure 17.2: Capturing region of adsorbed atoms of a one-sided model in front of 
a step (a) for a straight step and (b) for a curved step, (c) Capturing region for a 
symmetric model for a curved step. 


is shown in Fig.17.2. The area of the capturing region increases approximately by a 
factor 1 + x a /2p, and the velocity increases with the same factor. Thus the step with 
a curvature k= 1/p has a velocity higher than the straight one by 

&v A = v 0 ■ —ft. (17.5) 

The bump is accelerate compared to the straight part and is pushed further forward. 
The diffusion causes a destabilization of a step profile. 

Competing with this destabilization is the stabilization effect by the step stiff¬ 
ness. Since the equilibrium density increases at a curved part, the deposition flux to 
maintain equilibrium at a curved step increases as / eq (l + //lhK/k n T), and the driv¬ 
ing force or the supersaturation decreases correspondingly. The velocity of a bump 
decreases by 

6v e = -fe q ~;KX B Q 2 . (17.6) 

KbI 

Both destabilizing and stabilizing effects are, in the first order approximation, propor¬ 
tional to the curvature, as is apparent in Eq.(17.5) and (17.6). Since the diffusional 
instability increases with the velocity Vo, the instability wins eventually by increasing 
the deposition rate. The instability takes place when Sv A + Sv a = 0. The critical 
deposition rate / c is determined from the Eqs. (17.4-17.6) as 


h 


/eq 



2(Kh \ 

x<,kijT ) 


(17.7) 


Exercise: Explain that in the BCF model where K + = A'_ = oo, the step instability 
does not take place. 
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Answer: Consider a part of a step pushed forward with a positive radius p as 
shown in Fig. 17.2c. The velocity contribution from the lower terrace increases as 
v + = u 0 (l + x B /2p), but that from the upper terrace decreases as V- = v 0 (l - x s /2p). 
The total velocity » + + V- = 2vq remains independent of the step deformation in the 

stationary approximation. 

17.2 Structure of an unstable step 

In order to describe the profile of a destabilized step, one has to treat the problem 
quantitatively and analytically. First we treat the stability of a straight step in a 
linear analysis [8, 186], The step deformation y = £{x,t) is decomposed in Fourier 
modes. In the linear analysis there is no coupling among modes, and the consideration 
of a single mode with a wavenumber q is sufficient: 


(( x , t) = Vot + a q e u,t cos qx. 

(17.8) 

Here Vq is the velocity of a straight step moving in y direction. If the amplification 
rate ui q of the deformation is negative, the amplitude of the deformation diminishes 
and the straight step recovers. On the contrary, if u> q is positive, the mode amplitude 
increases and the straight step is unstable. 

The normal direction of this deformed step (17.8) is given by 

n = ( aC/ftML « (qa^ 1 sin qx, 1), 
y/l + W/dx)* 

(17.9) 

the curvature is 22 

K =li + (dC/dx)>)V>* q2a < eU ' tcos<lx ' 
and the growth velocity is obtained as 

(17.10) 

v — (0, d^/dt) = (0, v 0 + iv q a q e w,t cos qx). 

(17.11) 

Therefore the normal velocity is written up to the first order of a q as 


v n = (n • v) = vo + w,a,e“ ,i cos qx. 

(17.12) 


The adatom density c(x, y; t ) is also modified with a wavenumber q in x direction. 
For y oo, the modification should decay and thus the density is written as 

c(x,y;t ) = Coo + (c^ - c oo )e“ (v_, ' 0<)/l * + 6c q e u,t cosqxe~ A,( - y ~ V0t) (17.13) 


in front of the step. The first two terms are the density distribution around the 
straight step. Since the density satisfies the diffusion equation in the stationary 
approximation, 

v 2 c+^yf = 0, 


(17.14) 
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Figure 17.3: Dispersion relation of a step deformation with Schwoebel effect. 


the damping coefficient A, in y direction satisfies the relation -q 2 + A 2 — x~ 2 = 0 or 

A, - \!<? + K 2 - (17.15) 

From the local equilibrium condition at the step, Eq.(17.3), the density deviation Sc g 
is obtained as 

i _ .. ( c oo — Ceq Ceqfl‘2/7 2^1 One! 

6c » = “°’(~-isrV' (1716) 

where the curvature k is approximated in the linear approximation (17.10). From the 
continuity boundary condition 

_ , ^ dc „ /dc dc \ 

«» - - A ^n. + J - (17-17) 


one obtains the relation 

Hj 1 («o + w q a,c‘ ,,,< cosga;) « A [———e - ^"” 0 *)/* 1 - A q 6c q e w,t cosga:e _A,K_ ” 0<) | 

L A’ 8 J 

~ — (Coo — c eq ) (l — — c u,i coaqx\ — D B A q 6c„e u,t cosqx. (17.18) 

x s \ x a ) 

By comparing the zeroth order term in a,, one obtains the velocity n 0 of the straight 
step (17.4). From the first order of a q , the dispersion relation between u q and q is 
obtained as 


OJq 


—AA 




«o(A, ■ 


■ x s 1 ) - D s iq A, 


W ^_2 

"k B T q ’ 


CeqA/^ 2^1 

feflT 9 J 


(17.19) 


which is shown in Fig.17.3. The first term corresponds the diffusional destabilization 
or 6vd in (17.15), and the second term corresponds to the energetic stabilization or 
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6v a in (17.6) for long wavelength or small q (< x s ). The dispersion relation (17.19) 
can be further rearranged as 

w, = -fi q • v q (17.20) 


with the mobility 


and the restoring force 


M, = D e n 2 2 Ag^(> 0) 


(17.21) 


"> = k ' ~ (£ ~ 0 (‘ - 
“ ^ !+ Tnr(ir‘) ,,+(17 22) 

where the long wavelength limit, qx s < 1, is used in the last approximation. The ef¬ 
fective stiffness is the force constant for the step recovery modified by the diffusion 
destabilizing effect as 


Ei(X_ 

! \ fe q 


(17.22) 


A* = £ i- 


x a k B T 




(17.23) 


Here / c is the critical deposition rate defined in the previous subsection, (17.7). At the 
instability point / c , the effective stiffness vanishes, and the step looses the restoring 
force to the straight form. Equation (17.23) shows that the step stiffness depends on 
the deposition rate /. For small / /3 e n is large and the step is stiff, while for large 
/ it is soft. This variation of the step stiffness is observed in the experiment on Si 
[125, 126]. 

Near the instability point / c , we have a small parameter 

« = (17.24) 

Jc Je q 

For positive e, the dispersion relation for long wavelength has a positive maximum 
at q m of order e l/2 as shown in Fig. 17.3. The maximum value of w, is of order e 2 . 
This mean that near the instability point, the step modulation has a large spatial 
extension of order and relaxes slowly in a time of order e~ 2 . In the weakly 
nonlinear case the dynamics of an unstable step can be analyzed by rescaling the 
variables and extract slow dynamics [199, 13]. One scales space and time as follows; 

X = e l t 2 x/x a 

Y = y/x s 

T = eH/r, (17.25) 

and introduce a dimensionless field u and an interface position (, which are expanded 

in terms of e as 


- Q 2 (c - Coo) = «o + eui + e 2 «2 + 1 
= eH = (Hq + « 2 //i -I- 


(17.26) 
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The diffusion equation is described in these variables as 

eu xx + uyy — u = 0. 


(17.27) 


Here and later in this section, the subscript Ior7 means the derivative by them. 
The boundary conditions are written as 


u a 


V)(l + e) + e 3 H T = u Y - e 2 u x H x 
V c e 2 H xx 


-V c (l + e) 


2 [\+€ 3 (H X ) 2 } 3 / 2 ' 


(17.28) 

(17.29) 


where V (: = 2Qlc,., l ft/k !i r rx„ is the step velocity at the critical point / c . By comparing 
each order of e, one obtains the solutions and relations as 


O(e 0 ) : 

«o — Age Y , Aq = —V c 

(17.30) 

O(e') : 

t»i =A l cr Y , = 1 + H 0 {X,T) 

no 

(17.31) 

0(e 2 ): 

M 2 = A 2 e~ Y + -Ai'XxYe = H 0 + ^Ho + ^o.xx 

+ H, 

(17.32) 

0(e 3 ) : 

«3 = A 3 e~' + ^A 2% xxY e v + g^i.xxxxC^ 2 + Y)e~ v , 



^ = Wo + Ih 3 + H l + H 0 H , + W xx + H 2 . 

(17.33) 


To fulfill the boundary conditions (17.28) and (17.29) consistently at the order 
0(f 3 ), H 0 should satisfy the differential equation as 


1 dH 0 Id 2 H 0 3 d*H 0 1 (dH 0 
K dT ~ 2 dX 2 + 8 ax' 1 2 { dX 


(17.34) 


The linear part represents the dispersion relation (17.20-17.22) in this unit. The 
nonlinear term ( dH 0 /dX ) 2 is the lowest possible term compatible with the transla¬ 
tional symmetry of the system. There cannot be terms containing Ho such as Hg 
or HgdH 0 /dX since the evolution should be independent of the absolute position of 
the step itself. The equation (17.34) is called the Kuramoto-Sivashinsky equation 
[118, 177] and is known to have spatio-temporal chaos, as shown in Fig. 17.4. The 
front has many hills and valleys. Valleys shift randomly to the left or to the right, 
collide with each other and annihilate. A hill widely spread splits and a new valley is 
formed randomly. The same chaotic behavior is observed in the Monte Carlo simu¬ 
lation of the step advancing in the adatom diffusion field, as shown in Fig. 17.5 [170]. 
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h 



V 



(a) (b) 

Figure 17.4: Time evolution of the Kuramoto-Sivashinsky equation, showing spatio- 
temporal chaotic behavior, (a) The height h and (b) the gradient v ~ dh/dx. 



Figure 17.5: Monte Carlo simulation of the time evolution of an unstable step growing 
from the diffusing adsorbed atoms [170]. It shows the spatio-temporal chaos similar 
to Fig.17.4. 





Part IV 

Diffusion-Limited Growth: 

Pattern Formation 

So fax we considered the effect of surface structure and surface kinetics on the crystal 
growth. But there are other processes which controls the growth. Crystallization 
proceeds in the following sequences: 

1. Atoms to be crystallized are transported to the crystal surface (= chemical 
diffusion), 

2. they arc incorporated in the crystal at the surface (=surface kinetics), 

3. and the released latent heat should be transported away from the crystal surface 
(= heat conduction). 

If all these processes are fast enough, ideal growth laws can be realized. But in 
reality, the slowest process governs the growth rate as a whole, and the deviation 
from the ideal linear growth laws explained in part I is expected. For the crystal 
with an atomically flat interface, the second, surface kinetic process is the slowest 
one and it controls the growth. The growth law is found to differ from the ideal one, 
as explained in Part III. For rough surfaces, kinetics is fast, but still the growth is 
different from the ideal behavior when the transport processes (1) and (3) are slow. 
In this part we consider the crystal growth and its morphology, when the growth is 
controlled by chemical or heat diffusion processes. 

Diffusion field induces instability in the growing interface, as, for example, shown 
in the previous section 17. This diffusional instability causes variety of patterns in the 
growth shape of the crystal. Examples are the fractal structure of the diffusion-limited 
aggregation (Fig.1.3), dendrite in the melt or solution growth (Fig.1.2), lamellar 
structure in the eutectic growth (Fig.26.2). These topics are studied in the present 
part. 

18 Diffusion Equation 

Since the mass or energy is conserved, the material or heat transport follows a diffusion 
equation [122]. To consider the problem concretely, we consider the crystal growth 
from the undercooled melt where the heat conduction evacuates the produced latent 
heat from the crystal surface. Since the melt is cooled at a temperature T x below 
the melting temperature 7m, the Gibbs free energy of the crystal per volume Gs is 
lower than that of the melt Gl as (3.1) or 

AG = G l - G s = 


(18.1) 
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Figure 18.1: Temperature profile around the flat interface of the crystal growing in 
the supercooled melt. 

where L is the latent heat per volume. This AG drives the material to crystallize. By 
the crystallization, however, the latent heat L released at the interface heats up the 
crystal surface to a temperature T i} higher than the far field value T x . Therefore, the 
driving force of crystallization at the interface reduces to A G t = L(Tu — T;)/Tm. For 
many metals or some plastic materials, the interface is rough at the melting point and 
thus the growth rate is proportional to AG*, following the Wilson-Frenkel formula 
(3.5) or 

V n = K t (T m - T.) (18.2) 

with the kinetic coefficient Kt■ The remaining supercooling T< - T x drives the heat 
transport and prohibits the interface from being heated up by the latent heat. The 
heat transport in liquid is described by the heat conduction equation 

c p~ = fcV 2 T, (18.3) 

where C p is the specific heat per volume, k is the thermal conductivity. When the 
crystal grows with a normal velocity V n the latent heat LV n is produced per unit area 
in unit time. This heat should be transported away by the heat flow in the normal 
direction n as 

LV n = -fc(n • V)T = -kd„T. (18.4) 

Here we restrict ourselves to the one-sided model such that the heat is transported 
only in the liquid. Extension to the two-sided or symmetric models is straightforward. 
The driving force AG = L(Tu — T QO )/T M is partitioned into the kinetic part AG; 
and the transport part A G t = L(T, — T X )/T M . The interface temperature T, is so 
determined that the velocity determined from the kinetics (18.2) agrees with that 
from the transport (18.4). 

The fundamental equations (18.2-18.4) is derived for the melt growth where the 
heat conduction controls the growth. Similarly, in the solution chemical diffusion 
controls the growth of crystal. In order to generalize the situation and to take up 
the essential features of the problem, we scale variables in dimensionless form. The 
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dimensionless diffusion field is defined as 

»(x,f)= T(X L ^ r °° | ( 18 .5) 

and it follows the diffusion equation 

|jT = Dv2u (1.8.6) 

with the thermal diffusivity D = k/C p . Boundary conditions are the conservation law 


V n = —Dd n u, (18.7) 

and the Wilson-Frenkcl law at the interface 


V n = K'{A -m- dn ), (18.8) 

where K' — K t (L 2 /C p T m ) is the new kinetic coefficient, which will be denoted K 
hereafter. 

A = T Tldf < 18 ' 9 > 

is the dimensionless undercooling, which is the undercooling normalized by the tem¬ 
perature increase caused by the latent heat production. u { = (T, - T^/LC* 1 is the 
value of the diffusion field at the interface, 


/ = 2_Zk_ 

LL/C P 


(18.10) 


is the capillary length proportional to the surface stiffness 7 , and k is the curvature. 


For the rough interface with many steps and kinks, the kinetics is expected to be 
very fast. In the limit of K —* 00 , the local equilibrium is realized at the interface as 

Ui = A - dn, (18.11) 

instead of the Wilson-Frenkcl law (18.8). In this case, the crystal growth is totally 
governed by the diffusion. 


19 Flat Interface 


When the crystal grows steadily with a macroscopically flat interface in z direction, 
the growth velocity V is expected to depend on the undercooling A. In the coordinate 
frame comoving with the crystal ( x,y,z' ~ z — Vt), the diffusion equation (18.6) is 
transformed as 

1 du _ 2 2 du 

77^7 = V U + Tx~l' 

D at lr> az 1 


(19.1) 
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where 



(19.2) 


is the diffusion length. For a steadily growing crystal, the diffusion field does not vary 
in time in the moving frame, du/dt — 0, or 


2 2 du 

v “ + ;^ = 0 


(19.3) 


Taking the origin of z' coordinate on the flat interface, the diffusion field is solved 
as 

u(z') = Aexp (19.4) 

under the far field condition u(z' = oo) = 0. The diffusion length l D is the thickness 
of the diffusion layer, and characterizes the spatial variation of the diffusion field. 
From the conservation boundary condition (18.7) one obtains the relation 

o a 

V = -Dd z ,u = D— = AV. (19.5) 

‘D 

Thus, A = 1. From the kinetic boundary condition (18.8), the growth velocity is 
determined as 

V = K( A-l). (19.6) 

Since U{ = A = 1, the interface temperature in the conventional unit is 7’, = 
T x + L/C p by utilizing Eq.(18.5): The interface is heated up by the latent heat. 
Since T, should be colder than the melting temperature T M for the crystal to grow, 
A has to be larger than 1. Even though the melt is supercooled at A > 0, the crystal 
with a flat interface cannot grow steadily for A < 1. The latent heat released heats up 
the interface so high that the heat conduction cannot transport it quick enough with 
small supercooling. In the case of local equilibrium (K — oo), the time-dependent 
solution of the one dimensional problem is exactly obtained, and the growth velocity 
decreases as £ -1 / 2 . 


Exercise: We consider the crystal growth with a flat interface under the local equi¬ 
librium condition. When the interface is moving in z direction as z = ((t), the 
growth velocity is shown [99] to decreases as t~ l/2 by assuming the scaling form for 
the diffusion field as 

u(zJ)=u (m 


(19.7) 


(1) Show that the time-dependent diffusion equation (18.6) is transformed to the 
ordinary differential equation for u(w) with a new variable w = z/£(t) as 
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(2) If u depends only on w, then show that 


((f) = ViDPt 

(19.9a) 

u{w) = 2 Pe p P e- Pm2 dw, 

J W 

(19.9b) 

where P is defined by the relation 

A = \ArPe p erfc( %/P), 

(19.10) 

with the error function defined by 

erf c(\fP) = -^= Jj-6 x *dx. 

(19.11) 


Eq.(19.9a) shows that the interface velocity decreases to zero as ( oc t '/ 2 . 


Answer: 


(1) By changing variables, derivatives arc related as 


du(z/C(t)) 

dt 

du(z/C(t)) 

dz 


d 2 u{z/C(t)) 
dz 2 


dw du _ z( du 
dt dw ( 2 dw 
dw du _ 1 du 
dz dw ( dw 
1 cPu 
( 2 dw 2 ' 


(19.12) 

(19.13) 

(19.14) 


where ( = d(/dt. Inserting the relations (19.12-19.14) in the diffusion equation 
(18.6), one easily ends with (19.8). 


(2) By rearrangement, Eq.(19.8) can be written as 


CC _ cPu/dw 2 = 2p 
D wdu/dw 


(19.15) 


The first term depends only on time t and the middle term only on w, and thus 
should both be constant, which is set 2 P. From the first equation, CC = 2 DP, 
the interface advances in proportional to f 1 / 2 as in (19.9). Integrating the second 
equality, one gets 

In ^-= —Pw 2 + const (19.16) 

dw 


or 


OM = Ccr ,, m ; 
dw 


(19.17) 
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With the far field condition, u(oo) = 0, it is integrated as 


u(w ) = — C f 

J tu 


dw. 


(19.18) 


Prom the local equilibrium condition (18.11) at the interface, z = ((t) or w = 1, 
the integral constant C is determined as 


A = -C J™ e~ Pw *dw = ~^<J^erfc(\/P). 
The conservation law (18.7) is now written as 



By inserting (19.15) and (19.19) in (19.20) one gets (19.10). 


(19.19) 


(19.20) 


20 Spherical Crystal 


When the crystal grows in a spherical shape, it emits heat in all directions, and it may 
sustain steady growth. But it will be shown that the steady growth is not possible in 
this case either. We study the time evolution of the radius R and its radial velocity 
R = dR/dt, as shown in Fig.20.1. From the symmetry the diffusion field u is expected 
to depend only on the radial variable r in a spherical coordinate. Then the diffusion 
equation is written as 

u « 0. (20.1) 


1 du _ ( d 2 2 d \ 
D dt i^3r 2 + r dr ) 


In the last equality the stationary approximation is used, where the relaxation of 
diffusion field is very quick compared to the shape variation of the crystal. The 
stationary distribution of the diffusion field takes the form 


u(r) = A/r 


( 20 . 2 ) 



V= - 


d_R 

dt 


Figure 20.1: Crystal growth in spherical shape. 
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with an integration constant A. The far field condition u — 0 at r —» oo is satisfied in 
(20.2). The growth velocity V is determined by two ways, from Eq.(18.7) and (18.8), 
as 





(20.3) 


Here the relation of the curvature k = 2/7? with the radius of the sphere R is used. 
The integration constant A is then determined as 


7? 2 ( A - 2 d/R) 
R + D/K 


(20.4) 


and the growth velocity as 


dR = D{ A - 2 d/R) = A / _ pA 
dt R + D/K K~ l + R/D \ Rj 


(20.5) 


If the crystal radius R is smaller than the critical radius p c = 2d/A, sphere cannot 
grow. If R is greater than />,,, the sphere can grow. For small radius, p c < R <. D/K, 
the kinetics resistivity A' -1 against the growth is dominant, and the growth velocity 
A'A(1 - Pc/R) corresponds to that obtained phenomenologically in Eq.(7.2). When 
the sphere grows large with R > D/K, the diffusions! resistivity R/D becomes 
dominant and the growth rate varies as 


V 





( 20 . 6 ) 


For sufficiently large sphere as R » p c , the growth rate V = R is inversely propor¬ 
tional to the radius R, and by integration one gets 

A 2 (t) = ADt +A 2 (0). (20.7) 

Asymptotically as t —> oo, the radius increases in proportional to tK 2 , and the velocity 
decreases as V = R oc Thus the steady growth of sphere is impossible. 


21 Parabolic Crystal 


In the analysis so far, neither the planar nor the spherical crystals can grow steadily if 
the diffusion controls the growth. Is there any shape which allows the steady growth 
of a crystal? Ivantsov [89] has shown exactly that the needle crystal with a parabolic 
tip can grow steadily, if the crystal interface is in local equilibrium with infinitely 
fast kinetics and the surface tension is neglected, namely the crystal interface is the 
isotherm at the melting temperature 7 m. 

If the parabolic crystal grows steadily as 


z = Vt — 


x 2 + y 2 



2 R 


( 21 . 1 ) 
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Figure 21.1: Crystal growth in parabolic needle shape. The latent heat is emitted 
dominantly near the tip. 


as shown in Fig.21.1, its height increases in proportional to the elapsed time At, 
whereas the width increases only in proportional to Therefore, the growth 

takes place mainly at the tip region, and the latent heat is emitted at the tip. The 
total amount of the latent heat released from the tip of parabolic dendrite down to the 
tail with a radius r is Litr 2 V per time. The height of the dendrite is r 2 /2R and the 
surface area is proportional to r 3 /R, precisely calculated to be nr 3 /3R. Therefore the 
average heating per area is 3 LVR/r, which decreases for a long parabola at r -+ oo. 
Heat will not accumulate and a steady growth is allowed in this shape. 

We now treat the problem more quantitatively. When the crystal is growing 
steadily in z direction with a velocity V, one transforms as usual to the moving frame 
of reference as ( x , y,z' = z— Vt). For the present shape it is convenient to transform 
further to the parabolic coordinate (Fig.21.2) 


* = r-z', 

7) = r + z', 

6 - arctan(j//x) (21.2) 


with r = \Jx 2 + y' 1 + z n . (See Appendix A21). Then the diffusion equation in the 
steady state is described as 


1 {d du d 1 d 2 u 1 1 t du ,.du\ 

vTt \dri V dr) + + + = 


(21.3) 


where l D = 2 D/V is the diffusion length defined in (19.2). 

By denoting the interface as i] = ??,((, 0, t), the conservation boundary condition 
(18.7) is written as 


„ , r dr), n, + i di}. 


It , ( du 
•=~ h {’'%-• 


dr)i du 


Vi + £ drg du \ 

4 de 89)' 


(21.4) 
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Figure 21.2: Parabolic coordinate. 


and the local equilibrium condition (18.11) is written as 

u(r]i) = A, (21.5) 

since there is no capillary effect (d = 0). 

The paraboloid crystal of revolution (21.1) corresponds to the interfacial shape 
with a constant r? as = R. Due to the symmetry, the field u depends only on ij, 
and the diffusion equation reduces to 


_9 

9t) 



( 21 . 6 ) 


It is easily solved as 

u{ri) = A + C P Tf-'e-^dt). (21.7) 

J fl 

Here the local equilibrium boundary condition (21.5) is used. From the far field 
condition that u(r] —> oo) = 0, the integration constant C is determined as 


C = - 


A 

/k° r\~ x er r '^' a dt]' 


( 21 . 8 ) 


From the continuity boundary condition (21.4) with rji = R and di ],/= dru/dO = 
dr)i/dt — 0, one gets the relation 

R=-l 0 Rp. (21.9) 


Inserting Eq.(21.7) and (21.8) into the relation (21.9), one gets the Ivantsov relation 


r OO 

A = (fl/f D )e H/,D / r]- 1 e~’ ,/lD dr] = Pe p E l (P) (21.10) 

J ft 
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with the Peclet number P defined as 


R _RV 
h~ 2D' 


( 21 . 11 ) 


Here Ei(P) =/®i 'e X dx is the exponential integral function [1], In two dimen¬ 
sions, one obtains the steadily growing parabolic crystal 

(2Li2) 

with the Ivantsov relation 

A = \fitPe p erfc( y/P) (21.13) 

with the error function defined in (19.11) [83], (See Exercise). For a small supercool¬ 
ing A, Peclet number P is small and the Ivantsov relation is approximated as 


A 


P(-In P- 0.5772 •••) 
VttP 


for 3 dimensions 
for 2 dimensions, 


(21.14) 


and for high supercooling as A 
relation is approximated as 


1, the Peclet number P is large and the Ivantsov 


f 1-1/P 

\ 1 - 1/2P 


for 3 dimensions 
for 2 dimensions, 


(21.15) 


Horvay and Cahn [83] extended the Ivantsov solution to the elliptic paraboloid: 


z~Vt = ^(x 2 + a 2 y 2 ) + f ( 21 . 16 ) 

with the aspect ratio (ratio of x axis to y axis) a. The undercooling A and Peclet 
number P = VR/2D is related as 


1 p e p f°° e — [°° Pe w 

a \Jw[w 4- (a~ 2 - 1 ) P] h iJ(P + w){P + a 2 w) 


dw. 


(21.17) 


For a = 1 the shape becomes the three-dimensional paraboloid of revolution and 
the Ivantsov relation (21.10) is reproduced. For a = 0 the system reduces to the 
two-dimensional problem, and Eq.(21.17) reduces to Eq.(21.13). 

The Ivantsov solution indicates that the steady growth is possible in the form of a 
parabolic dendritic crystal. But the Ivantsov relation (21.10) shows that the growth 
velocity V and the tip radius R are not determined uniquely for a given undercooling 
A. It only determines the product VR ~ P. On the contrary, experiments show that 
for a given A, V and R are determined uniquely. There is some factor missing in the 
analysis by Ivantsov. 
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Exercise: Find the Ivantsov relation (21.13) for a two-dimensional dendrite, when 
the surface tension is absent. 

Answer: Take the 2D parabolic coordinate 


t = r-z', x = s/tj) 

r) = r + z\ z' = ±(77 - 0 


(21.18) 


with r = s/x 2 + z' 2 = (( + rj)/2, the diffusion equation in steady state is written as 


a 


du 






du\ l f du 
+ 7 - [Th 


, du 


d(j'hVd v 


o. 


(21.19) 


Boundary conditions at the surface tj = r), are expressed as the local equilibrium 
condition 

u(ifc) = A, (21.20) 

and the continuity condition 

. dij, T), + £ dr}, 


d£ + 2v dt 


. du di), du 


( 21 . 21 ) 


By taking the parabolic shape a s r]i = p, the diffusion field u depends only on T] from 
the symmetry and the diffusion equation becomes simple: 


d_ 

dt] 



i. _a«. 

+ T^Tr, ) ‘ a ' 


( 21 . 22 ) 


The solution which satisfies the far field condition u{t} -+ 
equilibrium condition ( 21 . 20 ) at t] = p is obtained as 


u(ri) = A 



Ip r] l t 2 e T, t ,D dri\ 

17“ 1 / 2 e'*/ ,D dTj j 


The continuity boundary condition, (21.21), is written as 


oo) = 0 and the local 


(21.23) 


, du A p 1/2 e p/lu 


(21.24) 


By changing the integration variable to x = \Jr]/l D and using the Peclet number 
P = p/l o, Eq.(21.24) is transformed to the desired result: 


A = 2 VTe p r c-^dx = VwPe p evk(VP) (21.25) 

Jy/T 


with the error function (19.11). 
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(a) (b) (c) 


Figure 22.1: Isotherm distribution in the supercooled liquid around the interface (a) 
for a flat interface, (b) for a deformed interface without capillarity, (c) and for an 
interface with Gibbs-Thomson effect included. 

22 Stability of a Flat Interface 

In the section 19 we discussed the steady growth of a flat interface and found that 
for A < 1, no steady state is possible. If the interface remains to be flat, the growth 
will slow down, and the growth velocity becomes zero even with a finite driving force. 
Another possible scenario of the time evolution of the flat interface is that it looses 
stability and deforms. We now study the latter possibility. 

In order to discuss the extreme case of the diffusion limited growth, we assume 
hereafter an infinitely fast kinetics at the interface, K -> oo, and the local equilibrium 
condition (18.11), or 

U{ = A — dK. (22.1) 

First we give a qualitative explanation on the stability of the interface. If the flat 
interface moves steadily with a velocity V, the diffusion field u varies from A to 0 in 
a length scale of diffusion length Iq = 2 D/V. The isotherm u — const lies parallel to 
the interface z' = z — Vt = 0 as schematically shown in Fig.22.la. When a part of 
the interface advances faster than the other by some fluctuation, what will happen to 
this pointed part? First we neglect the small effect of surface tension or capillarity. 
Then at the interface the diffusion field u, takes a constant value A, irrespective of 
its deformation. Near the pointed part the isotherm «,• = A pushes forward other 
isotherms, and the isotherm density becomes high here, as shown in Fig.22.lb. The 
high density of isotherm means the steep slope of Vu, and the heat flux —DV u 
increases and releases the latent heat quickly. The pointed part can grow faster than 
the undeformed part, and the deformation is enhanced: The flat interface becomes 
unstable. This instability is first studied analytically by Mullins and Sekerka, and is 
called Mullins-Sekerka instability [141, 122], 

So far we considered the interface instability when the crystal is growing in the 
supercooled melt and the heat is transported through the liquid. If the crystal is 
supercooled and the melt is hot, what will happen for the interface deformation? 
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Figure 22.2: Isotherm distribution when the crystal is supercooled, (a) A flat inter¬ 
face, and (b) a deformed interface without capillarity. 


Since the crystal is supercooled, the latent heat is released through the crystal. For a 
flat interface isotherms in the crystal are running parallel to the interface, as shown 
in Fig.22.2a. If a part of the crystal grows fast at some point by fluctuation, the 
separation between isotherms near the pointed part increases and the isotherm density 
decreases (Fig.22.2b). The slope of the isotherm decreases, and the latent heat is 
less effectively transported. Thus the growth rate there decreases, and the interface 
recovers the flat profile. In this case the interface is stable against fluctuation. There 
is an asymmetry and the instability occurs when the interface is propagating in the 
region where the transport is taking place. 

We now consider the effect of surface tension on the stability of a flat interface. 
The surface tension lowers the degree of supercooling at the pointed part with positive 
curvature (k > 0) as (22.1) and shown in Fig.22.1c. Since the driving force at the 
pointed part decreases, the growth rate there decreases. Thus the surface tension acts 
as a stabilizing factor for the flat interface. The total stability is determined by the 
competition between the diffusional destabilization and the energetic stabilization. 

We now describe the stability in the linear analysis. The interface is assumed to 
be deformed sinusoidally with a wavelength A = 2ir/q as 

( (x,y,t) = Vt + ft, exp (w,<) cos (qx). (22.2) 


The amplification rate w, determines the stability of the flat interface: When w, is 
negative, the interface is stable, whereas when w, is positive, it is unstable. Since 
the interface deformation influences the diffusion field u, it also has an additional 
variation to that of the flat interface in Eq.(19.4) as 


u(x,y,z-t) = A exp 


2 (z - Vt) 
fn 


+ 6u v exp(uj q t) cos(qx) exp [—A,(z — Vt)] . (22.3) 


Here A, should be positive, since the far field u(z —* oo) will not be affected by the 
interface deformation and the deviation should decay in z direction. By inserting 
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Figure 22.3: Dispersion relation of a sinusoidal deformation of a flat interface. 


(22.3) in the diffusion equation (18.6) and by considering the terms containing cos qx, 
one gets the relation 

/. i VA 

(22.4) 


6 J, 2A, 

D h ~ q 


2 + A 2 - 


From the local equilibrium condition (22.1) with the curvature k « -d 2 z[d: r 2 , A 
should be equal to A and the deformation amplitude Su q and a q are related as 


6u q = (^-A - d 0 q 2 ^ a q . 


(22.5) 


In the linear analysis the anisotropy in surface energy 7 gives no contribution and do 
denotes the isotropic part of the capillary length defined in (18.10). From the energy 
conservation or continuity equation (18.7), one gets the relation that A = 1 and 


= ~(B + (i- doq2 ) A<i - 

From Eq.(22.4) and (22.6) the damping in z direction is determined as 


15 


Id 

and the amplification rate as 
2 


A, = f-^W + Mti-^ + 


2d 0 , (%q 2 


Id 




( 22 . 6 ) 


(22.7) 


( 22 . 8 ) 


which is depicted in Fig.22.3. For slow growth, the diffusion length l D — 2D/V is 
large; Id/ do » 1. For a deformation with a wavelength less than l D or ql D > 1, the 
dispersion relation ( 22 . 8 ) is approximated as 


(22.9) 
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For wavenumbers q larger than the stability value 


q ‘~ V doin' 


( 22 . 10 ) 


the deformation damps down by the negative amplification rate u q . For the small 
wavenumbers q, however, the capillarity is not strong enough to stabilize the fiat 
interface, and u q is positive. In terms of the wavelength, the mode with a modulation 
wavelength longer than the stability length 





( 22 . 11 ) 


is amplified. The most unstable mode has the wavenumber 

q m = <J a /V3 


( 22 . 12 ) 


with the amplification rate 

= \vq m (22.13) 

as is shown in Fig.22.3. So far we considered the one-sided model. If the diffusion 
constant in the crystal is the same with that in the liqu id p hase, namely in the 
symmetric model, the stability length is given by A s = 2ny/d n ln- 


Exercise: Crystal in a spherical shape is growing in a static diffusion field V 2 u = 0 
with the local equilibrium boundary condition (18.11). Discuss its stability against 
the deformation as 

r(M) = R(t) + a lm (t)Yi m (0,<l>) (22.14) 

with Y/ m a spherical harmonics in spherical coordinates ( r,$, <j>) [140]. The Laplacian 
in spherical coordinate is expressed as 


V 2 


1 d 

i 

to 

r 2 dr 

\ d r . 

1 d 


) 

to 1 

1 

\ dr 


1 d f d\ 1 d 2 

r 2 sin 8 88 \ m 88 ) r 2 sin 8 d(f> 2 

+ 4>)i 

r 2 


(22.15) 


and spherical harmonics are eigenfunctions of the angular part A with an eigen¬ 
value — 1(1 + 1) as 

kY lm (6,<i>) = -l(l + Wm- (22.16) 

Also in the linear approximation the curvature k(6, 4>) can be represented as 
k{ 6, <j>) = | - + 2 )(l - 1 )a lm Y lm (8,4>). 


(22.17) 
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Answer: The diffusion field is modified from Eq.(20.2) as 

u(r,0,</>) = ~ + 6ui m (r,t)Y, m {6,4>). (22.18) 

Since the field u(r, 0,<t>) should satisfy the diffusion equation in the static approxima¬ 
tion V 2 n = 0 or 

d_^d6^_ l{l + l)6ujm = Qt (22.19) 

the deformation amplitude which decays far from the crystal is solved as 

fiutm = j.i+1 • ( 22 . 20 ) 


The local equilibrium boundary condition (18.11) at the interface (22.14) is explicitly 
written as 

Ui = aQ- %£r,m) + j^Y lm = A - d [| - ±(l + 2 )(l - l)a im n m ] (22.21) 
in the linear approximation, and two parameters A and B are determined as 

A 


-M) 


A- 


B = R 1 - 1 [A-(l + 2 )(l - 1 )d\ aim = R‘ 

The continuity equation(18.7) is expressed as 

Y n — R + CtlmYlm- 

From the zeroth order, one gets the growth velocity 


1(1 + l)d 


R 


aim- 


R\ R) R \ R! 


From the linear order in <Xf m , the modulation velocity is obtained as 


aim — %(l 1 ) 


A - ^(l 2 + 31 + 4) 
H 


0>lm — bJlm&lm • 


( 22 . 22 ) 

(22.23) 

(22.24) 

(22.25) 

(22.26) 


For small radius R when w< m < 0, the deformation amplitude ai m decays and the 
spherical shape is stable. For large radius when ui m > 0, the deformation grows and 
the spherical shape becomes unstable. The critical radius for the 1-th mode is given 
from the relation W| m = 0 as 

Rc.il) = + 31 + 4) = p f - 4 . (22.27) 

Since the l = 1 model is marginal as u lm — 0, the l = 2 mode is destabilized first at 
the radius R seven times of the critical nucleation radius p c = 2d/A. Thus a small 
sphere loses stability in the shape, and starts to deform. 



82 


Part IV. Diffusion-Limited Growth: Pattern Formation 



Figure 23.1: Fractal dendrite of Au deposited on Ru [86]. 

23 Fractal Dendrite 

If there is no surface tension (do = 0), the dispersion relation (22.8) reduces to 

£ - rkl, (23.i) 

U l n 

which is positive for all deformation wavenumbers q. The interface is always unstable, 
and it is most unstable for the largest |g| or the finest structure of deformation. From 
Eq.(22.27) the spherical crystal is also unstable for any size if there is no capillarity 
effect. Therefore, one expects that the very fine and irregular structure will be realized 
for the crystal growing in the diffusion field without the surface tension. This kind 
of structure is actually observed in the vapor deposition on a very cold substrate, as 
shown in Fig.l.3a [41] and Fig.23.1 [86], The adsorbed atom can diffuse randomly on 
the substrate, and its life time is very long on a cold substrate. In the meanwhile, 
the adatom collide with another adatom, coagulates and stop moving. Since the 
substrate is so cold that the once coagulated atoms never dissociate again to search 
for energetically more favorable position, the surface energy cannot play its role to 
stabilize the interface morphology. The irregular structure so grown is self-similar 
and called fractal, and is intensively studied recently [189, 10]. There are many 
other examples of fractal objects as aggregates grown by electrochemical deposition 
[130, 21] and bacteria colonies [132, 22]. We summarize the fractal theory relevant to 
crystal growth, and discuss fractal-to-compact crossover for the aggregate growing in 
a diffusion field of a finite density [191, 185]. 

23.1 Diffusion-limited Aggregation (DLA) 

The irregular aggregate is first studied theoretically in the computer [197]. An im¬ 
mobile crystal particle is placed at some point as a seed, and a randomly diffusing 
gas particle is released far from the seed. When the gas particle comes into contact 
with the seed particle, it freezes and is incorporated into the aggregate. Then, the 
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Figure 23.2: Diffusion-limited aggregation (DLA) grown in computer [150j. There 
are (a) 6 x 10 4 , (b) 6 x 10 5 and (c) 6 x 10 6 particles in each cluster. 


next gas particle is released far from the aggregate, and performs the random walk. 
After it makes contact with the aggregate, it freezes, and a new particle is released 
again. By iterating this procedure, the aggregate grows by incorporating atoms one 
by one. The grown aggregate is called the diffusion-limited aggregate (DLA), and 
its structure is shown in Fig.23.2. It is irregular, ramified and very open. The part 
looks similar to the whole by enlargement, and there is no characteristic length in 
the structure except the lower and upper cutoff lengths, the atomic and the system 
sizes. This self-similar object is called fractal. If there are N(r) atoms in a region 
with a radius r, the number of atoms in the radius br is N(br) = b D 'N(r). By taking 
b = r~ l , then 

N(r) = r D <N( 1) (23.2) 

and D( is called the fractal dimension. If the object fills the d dimensional space 
homogeneously with a finite density, Df = d and it is called compact. DLA is, on 
the contrary, fractal with Df = 1.71 in d — 2 dimensions, and with D f = 2.49 in 
d = 3 dimensions [134]. Since there are N ~ r D< number of particles in a radius r, 
the aggregate density is 

n(r) ~ p ~ r D, - d , (23.3) 

which becomes zero asymptotically for r -+ oo. It means that the aggregate has many 
open spaces. Zero asymptotic density is quite imaginable because the aggregate grows 
from a diffusion field of zero density: There is only a single gas particle during the 
whole growth simulation. Also a particle can walk forever until it comes to contact 
to the aggregate, and it means that the aggregate grows very slowly, V —* 0. The 
diffusion equation thus reduces to the Laplace equation V 2 u = 0. The diffusion 
length l d = 2 D/V is infinity, and the surrounding field as well as the aggregate have 
no characteristic length in the DLA growth. 
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Figure 23.3: Schematic behavior of an aggregation density n(r) growing in a diffusion 
field of a finite gas density n g . The diffusion length I D is the characteristic length 
scale of the crossover between the fractal and compact structure. 

23.2 Fractal-to-Compact Crossover of an Aggregate in a 
Finite Gas Density 

In a real experiment an aggregate grows from a gas with a finite density n 9 , whose 
time evolution follows the diffusion equation. The growth rate of the aggregate V 
should no vanish, and the relation between the velocity V and the gas density is 
of interest [185]. 

To realize the steady growth easily, we consider a unidirectional growth of an 
aggregate from a linear seed, similar to the “fractal forest” shown in Fig.l.3b. In 
front of the crystal, there are gas particles with an average density n g and they are 
making random walks on a square lattice to the nearest neighbor site. When a gas 
particle comes into contact to a seed or to an aggregate grown from it, it freezes 
and becomes a member of the aggregate. The time scale is so chosen that in a 
unit time every gas particle makes one diffusions! jump on average. By taking the 
lattice parameter as a unit of length, the diffusion constant D is 4. The aggregate 
front h(t) is defined as the height of newly crystallized particle, and the velocity V 
is obtained from the time variation of h(t). Since the gas particle is incorporated in 
the aggregate instantaneously and irreversibly, the density in front of the aggregate 
vanishes. It relaxes back to the given density n g within the distance of diffusion length 
Id = 2 D/V. Within this length, the gas density is low and the situation looks similar 
to the DLA growth: A gas particle sticks to the aggregate randomly and never melt 
back. Therefore, the grown aggregate has a fractal structure. But, since the aggregate 
grows from a gas with a finite density n g , the material conservation does not allow the 
aggregate density to diminish as n(r ) ~ r D, ~ d like the fractal object in the previous 
subsection does. For large distance r —* oo, the aggregation density should saturate 
to n 9 , as shown in Fig.23.3. A crossover of the aggregation density from the power 
law decay (23.3) to saturation takes place within a characteristic length Id, and the 





23. Fractal Dendrite 


85 



Figure 23.4: Irreversible and unidirectional solidification of aggregate from a gas with 
finite density: (a) n a = 0.08, and (b) n g = 0.1 [185]. 



Figure 23.5: Velocity of the aggregate front v versus the gas density n g , showing 
the power law dependence v ~ n g with u — 1 /(d - D { ) in d = 2 dimensions with 
Df = 1.71 [185]. 


relation 

- Ug (23.4) 

is expected. There should be a crossover from the fractal to compact structure at 
the length scale of order l D . This gives the relation between the aggregate growth 
velocity V and the gas density n g as 

V ~ J5 1 ~ n^ i ~ D ' ) . (23.5) 

In two dimensional system (d = 2), Monte Carlo simulation shows that the aggregate 
consists of many branches of irregular dendrites, and the separation between branches 
decreases as the gas density n g increases, as shown in Fig.23.4. This is compatible with 
the decrease of the characteristic length as n g increases. The growth velocity of the 
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aggregate V is found to increase in powers of the gas density n g with the exponent 
compatible with the fractal dimension Dt = 1.71 in two dimensions, as shown in 
Fig.23.5. The scaling (23.5) is recently confirmed in the deposition experiment of 
silver metal leaves from AgNC >3 solution [139]. 


24 Capillary Effect and Regular Dendrite 

Without the surface tension, the crystal profile is shown to take a random and irregu¬ 
lar form. In reality, most metals and some plastic crystals growing in a free and open 
space takes a regular dendritic form with a stable tip oriented in a special crystallo¬ 
graphic direction, as shown in Fig.24.1 [84], The tip of the dendrite is parabolic and 
the tip radius R and the growth velocity V satisfies the Ivantsov relation (21.10) with 
the undercooling A (Fig.24.2) [84]. In experiments the undercooling A determines 
R and V uniquely, contrary to the Ivantsov solution where an infinite degeneracy is 
expected. Capillarity may select a unique operating point among multitude of possi¬ 
bilities. It introduces a new characteristic length, the capillary length d = x /C p T^ / L' 2 
of (18.10). It is intrinsic, since it is determined only by the material parameters and 
every other length as p or l D is expected to be scaled with d. 

Since the main interest lies in the interface structure, it seems advantageous to 
integrate out the diffusion field and to derive the interface dynamics. Many local 
models for the evolution of the interface profile were proposed, and gave contributions 
to understand the problem of the velocity selection [40,17,18,100,101]. But the true 
interface dynamics is nonlocal even in quasi-stationary approximation, as explained in 
detail in Appendix A24.1. We summarize the results of recent studies on the nonlocal 
model. 

24.1 Velocity Selection by Microscopic Solvability 

Due to mathematical simplicity, the analysis is mainly done on two dimensional 
dendrite [124, 105, 33, 156]. The capillary length is defined in terms of the two- 
dimensional surface stiffness /?. The surface tension ft is assumed to have four-fold 
rotational symmetry reflecting the crystalline order of, such as, succinonitrile. 

P(0) = P o (l + €cos40). (24,1) 

For positive e, 6 — 0, ±7r/2, and n corresponds to the maximum of the surface 
tension, and to the corners in equilibrium shape (Fig.6.5). The surface stiffness is 
expressed as 

d 2 B 

P = d#Z = P°( 1 ~ 15 * C0B 4<? >, (24.2) 

and the capillary length (18.10) is then written as 


d = do(l — ecos 40) 


(24.3) 



GROWTH VELOCITY (cm/«J 


24. Capillary Effect and Regular Dendrite 



Figure 24.1: Dendrite of succinonitrile [84]. 



Figure 24.2: Velocity V versus the tip radius R for a given undercooling A [84] 
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with do being the magnitude of the capillarity length and e = 15e the strength of the 
anisotropy. Orientations 9 = 0, ±tt/2, n correspond to the minimum of the stiffness 
and the capillary length. 

By including the capillarity, the stability length appears, for example, A 9 = 
27ry'cWn/2 in the one-sided model. When the tip radius p is smaller than A s , the 
tip is unable to grow further due to the capillary effect, and the interface turns to 
be flat. When p is larger than A s , the tip is unstable to the deformation due to the 
Mullins-Sckerka instability. Therefore, the dendrite tip radius is expected to be of 
order A s . Then, one introduce a dimensionless parameter 

_ dpln _ do_ _ 2D dp _ d 0 v . 

° - p 2 “ pP ~ vp 2 ~ 2DP 2 ’ 1 ' ’ 

which is called the stability parameter. Since the capillary length do is small, this a 
is small, and one may try the expansion of the profile in powers of a. However, as 
the capillary effect is contained in the Gibbs-Thomson condition, it is coupled to the 
curvature k, the highest derivative of the interface profile as <jk -ad 2 z/dx 2 . In 
this case, normal perturbation fails and one has to use a singular perturbation [16]. 
There are many theoretical works on the dendrite theory which are summarized in 
many reviews [124, 105, 33, 156], 

In solving the nonlocal interface equation, boundary conditions has to be con¬ 
sidered appropriately. Far down the dendrite tail, the curvature is so small and the 
profile should approach to the Ivantsov parabola. By integrating the interface profile 
from the tail to the tip, one gets in general a finite slope (-) at the tip and the solution 
cannot be extended symmetrically to the other side of the tail. Slope at the tip can 
vanish only for special choice of the growth velocity v or the stability parameter o\ 
This is called the solvability condition [124, 105, 33, 156], and is sketched in a little, 
more detail in Appendix A24.2. Here! I summarize the main conclusions. 

1. Without surface anisotropy c = 0, there is no steady state for the dendritic 
growth. 

2. With an anisotropy c / 0, the symmetric dendrite grows steadily with its tip 
oriented in the direction of the stiffness minimum 9 = 0. 

3. The stability parameter a is found to depend only on an anisotropy parameter 
e, but is independent of the undercooling A. For small e, a is approximately 
proportional to c 7 /' 4 and thus vp 2 remains constant for various undercooling A: 

vp 2 = 2Dd n /a(e) oa e 7 ^ 4 . (24.5) 

The value of a{e) for the one-sided model is shown to be twice as large as 
that of the symmetric model [138]. By combining with the Ivantsov relation, 
vp = 2DP(A), the growth rate is obtained as v ~ P 2 . For d = 2 it reads as 
v « A 1 at small undercooling where the approximation (21.14) holds. 
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Figure 24.3: Simulation of the dynamical evolution of the dendritic crystal at (a) 
e = 0 (isotropic surface energy), and e = 0.10. The isotropic surface energy (a) shows 
the tip splitting, and the anisotropic one grows stably in (b) [165]. 



Figure 24.4: The stability parameter a dependent on the anisotropy parameter e at 
two different undercoolings A = 0.25 and A = 0.50 [165]. 


This result is confirmed by the numerical simulation by solving the diffusion equa¬ 
tion in a quasi-stationary approximation by the boundary element method and inte¬ 
grating the shape evolution [165], as explained in Appendix A24.3: For e — 0, the 
dendrite tip splits as it develops in time (Fig.24.3a), and for e > 0 the tip stable den¬ 
drite grows as shown in Fig.24.3b. The parameter a = d 0 v/2DP 2 is found to depend 
on anisotropy e but is independent of the supercooling A, as shown in Fig.24.4. 

In the three-dimensional experiment of succinonitrile dendrite growing in its melt 
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Figure 24.5: The stability parameter a at various undercooling A’s for succinonitrile 
[84], 

(Fig.24.1) [84], the relation VR 2 = const is found to be satisfied as shown in Fig.24.2 
and Fig.24.5. The strong anisotropy dependence of VR 2 is, however, not observed 
experimentally. There is some argument that anisotropy is not necessary in the 
mode selection [71], but at least it seems plausible that the dendrite grows in the 
direction of the stiffness minimum, since the minimum stiffness means less effective in 
suppressing deformation and the deformation grows fastest in that direction. If there 
is no preference iri orientation, the growing tip is very susceptible to the orientation 
fluctuation, and may split or form irregular dendrite. 

The chaotic behavior of the isotropic interface is already found in the Kuramoto- 
Sivashinsky equation discussed in section 17. There is another example of tip insta¬ 
bility in an isotropic system, a viscous finger problem [162, 152]: A viscous fluid is 
confined in a narrow gap between two glass plates of a Hele-Shaw cell [75]. When the 
fluid with low viscosity (like air) is pushed into the fluid with high viscosity (like wa¬ 
ter), the meniscus between two fluids is unstable and forms finger-like pattern, known 
as the viscous finger. The problem is formulated similar to the dendritic growth. One 
essential difference between the two problems is that the surface tension is isotropic 
in the viscous finger problem, whereas it is anisotropic in the crystal growth. The 
pattern realized in the radial Hele-Shaw cell shows a branched ramified structure as 
shown in Fig.24.6 [152]. By engraving a line groove in some portion on one glass 
plate of a Hele-Shaw cell, the tip stability is realized as shown in Fig,24,7a [131]. 
The fluid dendrite thus obtained satisfies the relation vp 2 — constant as shown in 
Fig.24.7b, until the tip radius p becomes too small. This experiment clearly shows 
the importance of the anisotropy in stabilizing the dendrite tip [19], 

So far, theories and simulations dealt with two-dimensional(2D) dendrite growth. 
Recently, solvability condition is successfully applied to the three-dimensional (3D) 
dendrite [14, 36], The surface free energy 7 (0,0) now has an orientation dependence 
on two Euler angles, 6 and <j> of the normal vector to the surface. In the theory, a 
cubic anisotropy 

7 (0, <f>) = 1 + e(4 cos 4 6 + 3 sin 4 9 + sin 4 8 cos 40) (24.6) 
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Figure 24.6: Viscous finger pattern in a radial Hele-Shaw cell [152]. 



0-00! 001 G 1 T 10 

TIP SPEED V (cm/s) 


Figure 24.7: (a) Viscous finger with a groove line cut on one glass plate. Along 
the groove grows a tip-stable parabolic, dendrite, whereas other tips are unstable to 
splitting, (b) Log-log plot of the tip radius of the fluid dendrite versus tip speed [131]. 
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is assumed. The selection problem is solved in two stages: Near the dendrite tip, the 
shape is perturbed as 

z = + A 4 r 4 cos40-|- (24.7) 

with r 2 = x 2 + y 2 . The growth rate V and the tip radius R is found to satisfy the 
same relations as in the 2D case: the Ivantsov relation 

VR = 2DP( (24.8) 

for small supercooling A, and the solvability scaling relation 

VR 2 = 2Dd,~ e -7/< (24 .9) 

o(e) 


Therefore, the anisotropy dependence of the velocity V and the tip radius R is ex¬ 
pected to be the same with the 2D case [14]. 

The analysis near the tip cannot be extended to the dendrite tail straightforwardly, 
since the correction A 4 r 4 cos 40 grows faster than the unperturbed Ivantsov paraboloid 
-r 2 /2R [36], But the correction form indicates that four fins extend in the directions, 
0 = 0, 7r/2, 7 t, 3tt/ 2. When one looks far down the dendrite, the crystal interface is 
almost parallel to z axis, and the diffusion field hardly varies in z direction; d 2 ujdz 2 « 
0. Then the evolution of the dendrite fin in the tail region is essentially controlled by 
the two dimensional diffusion [36]. The four fins that have started to grow near the 
tip develop into two-dimensional parabolic dendrites in xy cross section. For example, 
down at a height z(< 0) from the tip, a dendrite growing in x direction takes a form 

* = 3 -r (24.10) 


where p and v are the selected values of radius and velocity of the 2D dendrite. Here 
the time t is replaced by the height J^:| divided by the 3D growth velocity V, since 
the profile at a height z develops to that at a height z — Vdt after a time dt in the 
steady state 

Between the two extreme asymptotics, Eqs. (24.7) and (24.10), an intermediate 
asymptotics is found [36] until the arm length exceeds the 2D diffusion length. In 
this intermediate region four fins affect their growth mutually. At a height |z| from 
the tip, the total cross-sectional area of fins should be that of the Ivantsov paraboloid 
for the steady growth. If the fin protrudes a distance x with a width y, then the area 
is about xy ~ |^|. Even though x and y depends on the time \z\, the dimensionless 
parameter of the system 


<72 


vp 2 dx . 

2 Dd a d\z\ 



(24.11) 


may remain constant [5]. By assuming the profile x ~ \z\ a and y ~ \z\ l “, a-i is 
calculated as <t 2 ~ |^| 3 ' 5 “. If <x 2 remains constant, then a = 3/5, and the arm length 
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Figure 24.8: (a) Xenon dendrite and (b) tip contour compared with a power law fit 
to parabola and z ~ x 1 - 67 [26]. 

x grows as x ~ |^| 3 / 5 and the arm width y increases as y ~ |z| 2 / 5 . This behavior 
is originally found in the two-dimensional Hele-Shaw flow problem numerically and 
analytically [5], and later applied to the dendritic growth [36]. The profile of the fin 
in xz cross-section is written a s z ~ x b ^ 3 in this intermediate region. The profile 
is recently observed in the experiment of Xenon dendrite as shown in Fig.24.8 [26]. 
When the arm length exceeds the 2D diffusion length, the final aymptotics sets in 
and the arm length grows with a constant velocity. For the details, see the references 
(14, 36], 

24.2 Tip Stability and Sidebranches 

With an anisotropy e in the surface tension, the parabolic crystal with a tip radius R is 
found to grow steadily with a velocity V in the supercooled melt with the supercooling 
A. V and R are uniquely determined from both the Ivantsov relation VR= 2DP(A) 
and the solvability condition VR 2 = 2 Dd 0 o~ l (e). However, the dendrite is growing 
in the environment with fluctuations, for instance, created by thermal noise or by 
hydrodynamic convections. Is the tip stable against these fluctuations? Another 
problem is the origin of the sidebranches. According to the solvability theory, the tip 
of the selected dendrite is almost parabolic without sidebranches. How, then, are the 
sidebranches created in experiments? 

If the crystallizing front is flat, Mullis-Sekerka instability takes place and the 
front is unstable. It is most unstable against the fluctuation with the wavenumber 
q m = with the exponential amplification of the amplitude with the rate 

u m = 2vq m /3, as explained in the section 22. But actually the dendrite tip is curved 
parabolic as 

z =ax) = vt-^ 




(24.12) 






94 


Part IV. Diffusion-Limited Growth: Pattern Formation 



Figure 24.9: Trajectory z = £(x) of the deformation on the interface, which is alway 
perpendicular to the interface profile, z = £(x). 


in two dimensions. Since the interface deformation grows in the normal direction 
n — ( n x ,n z ) oc (—d(/dx, 1) as shown in Fig.24.9, the node of deformation follows the 
trajectory z = £(x) 

df 


or, by integration, 


/X') _p 
dx n x \dx) x 


£(x) = pin a- + C. 


(24.13) 


(24.14) 


When the initial deformation at t = 0 is given at position (xo, -x\/2p) on the 
interface, the integration constant C is determined as C = — x\/2p — plnxo- At time 
t, the deformation has been propagated at the site (x, ((x)) determined by crossing 
point of C(x) and £(x) as 


£(x) = vt - 


2 P 


Pin--/. 
x 0 2 p 


(24.15) 


Asymptotically for large t, x increases as y/Zpvt, but its height z = £(x) increases 
slowly as p/2 -In t. The vertical separation between the dendrite tip at (0, vt) and the 
deformed position at « ( x/2pvt, p/2 lu t.) increases as vt. Though the deformation 
amplifies with the rate o\„, it is convected down with a velocity v and the tip remains 
stable. This is called the convective stability of the dendrite tip. Down convectcd 
noise is amplified to form sidebranchcs. If this scenario is correct and sidebranches 
are originated from the random noise at the tip, sidebranches on the different side 
of the dendrite arc expected to be uncorrelated. By measuring the correlation the 
noise origin of sidebranches is confirmed [54]. The external oscillatory flow imposes 
systematic fluctuation at the dendrite tip, and the synchronized formation of side- 
branches is observed [28]. Of course, the tip stability depends on the amplitude of the 
initial noise. If the initial noise is too strong, the deformation is strongly amplified. 
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Figure 24.10: The sidebranch periodicity A near the dendrite tip is proportional to the 
wavelength of the most unstable mode A m ~ p\fo, irrespective of the undercooling A 
and the anisotropy e [165]. 


When the deformation amplitude becomes larger than the tip radius p before the 
deformation is convected down the length of order p, the dendrite tip splits [37]. In 
fact, the fractal structure is realized for NH 4 C1 dendrite growing in a Hele-Shaw cell 
when a bottom glass is cut rough randomly [82]. 

We now study the sidebranch formation more quantitatively. The convected de¬ 
formation amplifies and becomes a sidebranch with a periodicity A about that of the 
most unstable mode, 


= 27r^/|pv/o- (24.16) 

In the numerical simulation, the wavelength A of the sidebranch is defined by dividing 
the growth rate v by the number of sidebranches produced per unit time. The ratio of 
the sidebranch periodicity A to p^fcr was fount to be independent of the supercooling 
A and the anisotropy e as shown in Fig.24.10 [165]. 

As the perturbation slides down along the side of the dendrite z = —x 2 /2p, the 
normal velocity v n decreases as v„ = —v(d(/dx)~ 1 ~ v vWR. as shown in Fig.21.1. 

Then the local diffusion length increases as ln(z) = 2D/v n (z) ~ lv\J\z\/p as well as 
the most unstable wavelength 

A m (*) « A m (0) '* ~ lu(z ) l/2 . (24.17) 

The periodicity is expected to increase as the sidebranching is convected down the 
shaft of the dendrite. Along with the variation of the periodicity A m (z), the amplifi¬ 
cation rate u> m (z) = 2v n (z)q m (z)/3 of the most unstable mode varies while the 






Figure 24.11: (a-h) The time sequence of the disturbance amplification near the tip of 
the dendrite to the sidebranch. (i) The maximum curvature sq versus the arc length 
s [157]. K| is proportional to the noise amplification. 


perturbation is convcctcd down. 


u> m (z) 


2 v n {z)q,„{z) 
3 


v 

A m (0) 



(24.18) 


Therefore, when the perturbation is convected down to a height \z\ during a time 
t = \z\/v, the deformation is totally amplified by a factor A; 


' exp 


f u) m (z)dt 
Jo 


exp 


0 3 /' 1 


vl/4 


LA m (0) 


exp 


(2/3) 


1/4 


- 1/2 


27T 


kl 


1/4' 


. (24.19) 


Precise expressions of A arc obtained for the two-dimensional dendrites [155, 11], and 
for the three-dimensional axisymmctric dendrite [123]. If the initial perturbation is 
the thermal noise at the dendrite tip, it is found to be too small to explain the observed 
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position of the first sidebranch. On applying a noise near the tip by a localized heat 
pulse, the noise amplification to the sidebranch (24.19) is studied experimentally 
[157], The maximum of the curvature on the sidebranch is found to increase as 
k t ~ exp(as 1/4 ) during the noise is convected down the dendrite till the arc length a, 
as shown in Fig.24.11. 

As explained in the previous subsection, there is the latest theory of the three- 
dimensional non-axisymmetric dendrite with four fins formed a little down the tip 
[36]. Since the profile of the fin z = ((x) ~ r ,/3 has a slope as d(/dx ~ z 1 ^ which 
is milder than that of the parabolic dendrite dQv/dx ~ z 1 / 2 , the convective effect is 
small and the dendrite is more susceptible to the noise. At the height \z\, the most 
unstable mode has the short wavelength A m (z) ~ z l/r \ the large amplification rate 
w m (r) ~ z~ 3 / r \ and the large total amplification factor A ~ exp jcv^ 2 /^ 5 ]. The 
more precise formula is given by Brener and Temkin [38]. The amplification factor A 
increases steeply along the fin compared to the parabolic case, (24.19). The first side 
branch is in fact observed on the fin of the Xenon dendrite at the position expected 
from the thermal initial noise [26]. 

24.3 Dendrite in a Channel 

For the case of a free dendrite, the anisotropy of the surface energy is shown to 
stabilize the dendrite tip. When the crystal grows in a channel with non permeable 
walls, interaction with the channel walls through the diffusion field allows stable 
stationary patterns even for a system without anisotropy [104], The wall provides the 
effective anisotropy in the system and the steady growth of a symmetric finger-shaped 
crystal is possible if the dimensionless supercooling is large, A > i [31]. Due to the 
global conservation of the energy, the width w of a symmetric finger in a channel 
of width A should be equal to AA; w = AA. If A is small, the symmetric finger 
is far from the wall and it cannot maintain the stable form. Therefore, A should 
be larger than 1/2. For a fixed supercooling A, the growth velocity v varies as a 
function of the channel width A as shown in Fig.24.12. The result is obtained by 



Figure 24.12: Growth velocity v versus the channel width A [31, 35], 
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Figure 24.13: Asymmetric finger along the wall [35]. 

the extension of the numerical simulation described in Appendix A24.3. Velocity v 
has a maximum, and it decreases for a large channel width A. For an infinitely wide 
system, as we discussed section 24.1, the isotropic dendrite cannot grow steadily with 
a finite velocity. When the width becomes too wide, the finger splits to have shorter 
period. Solvability theory by Brener, Geilikman and Tcmkin [31] found the maximum 
velocity v ~ (D/do)(A — |) 7/2 at the channel width A ~ do{& — It is natural 

to assume that the state with a maximum velocity is selected and realized. 

In the simulation, however, the asymmetric dendrite is observed which grows along 
the wall, as shown in Fig.24.13 [35]. In the simulation a mirror boundary condition is 
used: u(x,y) = u(—x,y) — u(X — x,y) for the diffusion field and ((x,y) = ((-x,y) - 
((X — x,y) for the interfacial profile with 0 < x < A/2. This asymmetric dendrite 
along the wall thus really means a double-finger structure with its mirror image at 
negative x. This structure is called “doublon” for short [87]. In a time-dependent 
simulation model [87], the doublon is also found in the system with periodic boundary 
conditions. If two fingers are growing side by side, one generally anticipate that the 
one which steps little ahead by fluctuation wins the competition against the other 
by gaining more diffusion field supply. In the present situation, when one wins, the 
periodicity A increases and the growth rate v decreases according to Fig.24.12. The 
one ahead is caught up by the nearby finger. Thus the doublon can survive as a stable 
profile even at A < |. Recent solvability analysis shows that the doublon grows with 
the velocity v proportional to v ~ (D/do)A 9 without an anisotropy [15]. 

24.4 Morphology Diagram 

From the findings in the previous subsections, an isotropic crystal (if it exists) can 
grow radially outward with a doublon structure in an open space. Since the separation 
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Figure 24.14: (a) Equilibrium shape of a crystal coexisting with lattice gas atoms in 
a closed system. Simulation results (symbols) for different A//, are compared with the 
theoretical results drawn by curves. Time evolution of a crystal in an open system 
are depicted at different driving forces; (b) dendritic shape at Ap/fc/jT — 10, and (c) 
irregular structure with tip splittings at A r/UbT — 12 [166]. 

between two double-fingers widens as they grow, the space should be filled by new 
fingers developed from sidebranches. The structure as a whole should have a spherical 
and convex envelope, and is called the compact seaweed (CS) structure [34] or dense 
branched morphology(DBM) [20]. With an anisotropy, dendrites grow in directions 
of stiffness minimum. When the dendrite tips are separated radially, the secondary 
arms from the sidebranches fill the space with a concave envelope. This is called the 
compact dendrite (CD) [34]. 

The morphological transition from CD to CS pattern is first obtained in the Monte 
Carlo simulation of a crystal growth from a lattice gas [166], as is shown in Fig.24.14. 
The growth is similar to the DLA growth from a finite density gas, but with an 
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Figure 24.15: Morphological transition from (a) dendrite, (b) intermediate, and (c) 
dense branched morphology (DBM) of a hcxatic liquid crystal [151]. Undercoolings 
are (a) A = 0.24, (b) 0.45 and (c) 0.63. 

interfacial energy included. The interaction is the same with the Ising model (6.33) 
such that a broken bond from a crystal atom costs an energy J. When the gas atom 
diffuses and touches the crystal interface, it tries crystallization. If it crystallizes, there 
is a chemical potential gain -A/< but has to pay an energy cost (z/2 — n)J if the 
crystallized site has n nearest neighbor crystal atoms among z coordination number. 
These energy changes should be taken into account in the Boltzmann weight as is 
described in the subsection 9.1. One also has to consider the melting process from 
the crystal surface to satisfy the detailed balance. 

In a closed system where the total number of gas and crystal atoms is fixed, an 
equilibrium shape is realized as is shown in Fig.24.14a. In an open system, the gas 
atom is fed from a particle reservoir with a fixed density at a distance far from the 
growing crystal. Then the crystal grows steadily. For small Ap as Ap/fc s T = 10, 
the crystal grows in diagonal [11] direction in a regular dendritic form, as shown 
in Fig.24.14b. At a larger Aas A h/IcbT = 12, the crystal grows in an irregular 
form with concave envelope, as shown in Fig.24.14c. Thus simulation clearly shows 
a crossover in the growth morphology. There are many similar simulation works on 
the morphological transitions [175, 176]. There is also an experimental observation of 
the morphological crossover in the growth of columnar hexagonal crystal, as shown 
in Fig.24.15 [151], 

Analytical studies are also performed on the dynamical selection of morphology 
[34]. It is natural to assume a maximum velocity criterion that the pattern with 
maximum velocity is selected. With this hypothesis, Brener, Tcmkin and Miiller- 
Krumbhaar have derived the morphological phase diagram in the phase space of the 
supercooling A and the anisotropy e , as shown in Fig.24.16 [34], Since they didn’t 
considered doublons in an isotropic or weakly anisotropic region, CS structure can 
appear only for A > 1/2. By considering the double-finger structure, the crossover 
between the CS and CD structure is expected to take place at A ~ e 7 / 20 . 

At small undercooling A and small anisotropy e, the noise is expected to be 
important and induce tip splitting. When the tip is destroyed, the structure becomes 
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Figure 24.16: Morphological phase diagram in the phase space of the undercooling 
A and the anisotropy strength e. CS: compact seaweed, equivalent to DBM, CD: 
compact dendrite, FS: fractal seaweed, and FD: fractal dendrite [34], 

fractal for a short length scales, as is described in the subsection 23.2. Those struc¬ 
tures are called fractal seaweed (FS) and fractal dendrite (FD) in their morphology 
diagram, Fig.24.16. For the details, refer the original paper [34]. 

25 Unidirectional Solidification from the Solution 

An example in which the growth pattern is controlled by the material diffusion is 
the alloy crystal growth from solution. Solution of a binary alloy is encapsulate in 
a thin Hele-Shaw cell, a thin rectangular parallelepiped cell (Fig.25.1). One end of 
the cell is kept hot and the other end is kept cold to impose a temperature gradient 



Figure 25.1: Unidirectional crystal growth from solution in a Hele-Shaw cell. Solution 
is placed in a narrow gap between two glass plates, heated at one end and cooled at 
the other end. The cell is pulled steadily in the cool region, facilitating the steady 
growth of the crystal. 
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II 


kc L (T) 

Figure 25.2: Equilibrium phase diagram of a binary alloy. The liquidus temperature 
T l decreases from the melting temperature Tm of the pure A material by increasing 
the concentration c. of B atoms. 

over the sample. When the cell is pulled in the cold side, the solution is crystallized 
unidirectionally. The crystal-liquid interface lies normal to the pulling direction. 

Since the liquid is hot and the crystal is cold, the latent heat is released in the cold 
crystal. In this case, the heat diffusion will not induce instability of the flat interface 
as explained in the section 22. But there is a slow material transport which controls 
the crystal growth. The solubilities in the liquid and crystal are different, as is shown 
in the phase diagram Fig.25.2. The component less soluble in the crystal is expelled 
in the liquid and should be transported far away down by the diffusion. The other 
component, more soluble in the crystal should be supplied far away from the liquid 
by the diffusion. This material diffusion causes the similar interface instability as 
the heat transport did in the melt growth, and leads pattern formation with periodic 
structure. To summarize, in the unidirectional growth heat conduction stabilizes the 
interface and the chemical diffusion causes the instability. 

Since the heat conduction is faster than the material diffusion, and the heat con¬ 
ductivity in the crystal and liquid are almost the same, one often uses the approxi¬ 
mation that the thermal gradient Gt is the same in both phases. When the sample 
is sandwiched in a thin cell, the heat conduction is mediated by the cell wall, and 
the thermal conductivities in both crystal and liquid phases arc the same. Then the 
temperature distribution is simply described as 

T(z) = T 0 + G t z. (25.1) 


25.1 Fundamental Equations 

We consider the phase diagram of AB alloy as shown in Fig.25.2. By including B 
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Figure 25.3: Concentration distribution drawn in a heavy line at the vertical axis 
representing the temperature T and the height z simultaneously from the relation 
T = T 0 + GtZ. Here T a =T M - mc^/k. 


atoms with a concentration c, the melting temperature decreases from that of the 
pure A atoms, T M . At a temperature T below T M , the crystal phase is stable when 
the concentration c is less than cs(T), and the liquid phase is stable when c is more 
than c l (T). Between these two concentration, c s (T) < c < ci.(T), the crystal and 
liquid phases coexists. Phase separation between the crystal and the liquid phases 
takes place, and a crystal with a concentration c s (T) coexists with liquid with a 
concentration Ci(T). The liquidus line is well approximated by a line 

T l = T m - toc l (T) (25.2) 


with a slope m. The partition ( or segregation) coefficient k is defined as the concen¬ 
tration ratio between the coexisting two phases as 

,_ c s(T) 


Ci(T)’ 

and the miscibility gap is defined as 

Ac = c l (T) - c s (T). 


(25.3) 


(25.4) 


By pulling the two-dimensional cell with the solution of concentration c 00 to the 
low temperature side with a constant velocity V, the steady crystallization takes 
place. From the global conservation of material, the concentration of the crystal alloy 
should also be c x , as denoted in the equilibrium phase diagram, Fig.25.3. First we 
consider the case when the surface is flat at z = 0. Since the material diffusion in the 
crystal is negligible, the concentration in the crystal should be Coo = c°. Assuming 
that the local equilibrium is satisfied at the interface, the liquid concentration at the 
interface is c\ = c x /k and an interface temperature is 


(25.5) 
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A concentration difference is created in the liquid at the interface c® and far from 
there Coo. The concentration variation follows the diffusion equation 

f ^g + D^c, (25.6) 

with D ( being the chemical diffusion constant, and the first term on the r.h.s. is due 
to the fact that the laboratory frame is comoving with the crystal with the velocity 
V. 

When the interface deforms as 


* = COM). (25.7) 

the temperature at the interface varies as 

T, = T 0 + GtC (25.8) 


Since the Helc-Shaw cell is thin, we can neglect variation of physical quantities in the 
directions of the thickness, {/-direction. The equilibrium concentration of liquid 
at the interface temperature 7) is determined from Eq.(25.2) with the Gibbs-Thomson 
curvature effect included as 

T. = T M (l - — k) - mc. w , v (25.9) 


with L being the latent heat per volume. By assuming a local equilibrium at the 
interface, cr,(-r, () = cj„ ( , q , the liquid concentration at the interface is determined from 
Eq.(25.8-25.9) as 


MM) 


Coo qflTM G t 

-— -K - C 

k mL m 


(25.10) 


From the local equilibrium assumption, the concentration in the crystal cs(x,£) is 
given by cs(.r,C) = kc,i,(x,0 with the equilibrium segregation coefficient k. On crys¬ 
tallization, excess mass V„(<h.(M) - cs(.r, ()) is expelled from the crystal per unit 
time, and this excess mass has to be transported by chemical diffusion flux —D c dc/dn. 
The mass conservation reads as 


V n {c\. - c s ) = — D<. 


dc. 

dll' 


(25.11) 


We now introduce a dimensionless diffusion field in the liquid by 


u(x,z,t) 


C Coo 

Ac 


C Coo 
(k~ l - l)Coo 


(25.12) 


which satisfies the far field condition u(x, z —* oo) - 0. For a flat interface (£ = k — 
0), u is 1 at the crystal-liquid boundary. The diffusion equation reduces to 
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with the diffusion length l D = 2 D c /V. When D r . is very large, the diffusion field u 
relaxes to steady state du/dt = 0 quickly during the slow growth of a crystal. The 
field u follows the steady state equation 

— —+ V 2 « = 0. (25.14) 

lp oz 

The local equilibrium condition (25.11) is written in this dimensionless form as 


Ui y {. = 1 - d« - —. 

h 


(25.15) 


A chemical capillary length 

^ _ 7HT m 
mLAc 

characterizes the surface effect, and a thermal length 

mAc 
T ~ ~G^ 


(25.16) 


(25.17) 


characterizes the scale of the temperature variation. Since the concentration of the 
solid at the interface is 

m.-,s = - !). (25.18) 

the material conservation (25.11) is written as 


V n [k + (1 - fc)«i,Ll 


-D ^ 
"dn 


(25.19) 


Here V n is the normal growth velocity given by 



(25.20) 


25.2 Stability of a Flat Interface 

When the flat interface is growing steadily in 2 direction with a velocity V, the 
diffusion equation has a solution 

u 0 (z) = e~ 2 ^, (25.21) 

which satisfies the boundary condition u 0 (z = oo) = 0 and u 0 (z = 0) = 1. When the 
interface deforms to 

2 = ((x , t) — a q e ljl,t cos qx, (25.22) 

the diffusion field also deforms as 

u(x,z,t) = u 0 {z) + Aie"' 1 cosqxe~ A,z . 


(25.23) 
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Figure 25.4: Dispersion relation of the interface modulation during the unidirectional 
solidification. 


Here Ai is proportional to a,, and the damping rate in 2 direction, A,, should satisfy 
the diffusion equation in quasi-stationary approximation (25.14) as 

or _-- 

A, = Ip + V< f + l n 2 - (25.24) 

From the local equilibrium condition (25.15) A\ is determined as 


A ' = (h~h~ dq2 ) (25,25) 

By inserting these results in the material conservation (25.19), one obtains the dis¬ 
persion relation 

| = (<» + V?Ti?) (£ - i - *f) - ©’+£u - *> (r T +O ■ < 25 - 26 > 

which is depicted in Fig.25.4. For slow pulling velocity V, u> q is negative and the 
deformation decays in time. The flat interface is stable. When the pulling rate V is 
larger than the critical value 14, the maximum value of tj q becomes positive, and the 
flat interface becomes unstable against those deformations with wavenumbers with 
u> q > 0. At the critical velocity V c , the maximum of ui q becomes zero: ui q = du> q /dq = 
0. (See Fig.25.4.) Since 14 is small, the diffusion length l D — 2D jV is large such that 
the relation li j, / T d holds, and the wavenumber for the most unstable mode is 
expected as large as qlo > 1. Then the dispersion is approximated as 
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<?a 



Figure 25.5: Neutral curve of the unidirectional solidification. 


The critical velocity V c is obtained as 

K _ _2_ _ J_ J_ (4>fd\ 1/3 _ 3 /2k?d\ 1/3 _1_ 

D c f d,c h ^d,c \ Id,c ) h 2 \ It / h' 
and the critical wavenumber is 



(25.28) 


(25.29) 


For V > V c , u q is positive for some regions of wavenumbers, and the flat interface is 
unstable for these sinusoidal modes with wavevectors within this region. The locus 
of cj q (q, V) = 0 represents the neutral curve q = Qn(V) as shown in Fig.25.5. 

The neutral curve is shown to be closed at the upper critical velocity 14- When the 
pulling velocity is too fast, the diffusion length becomes of the order of the capillary 
length, It Id ~ d, and the interface is stabilized by the capillarity. This is called 
the absolute stability. For ql D < 1, the dispersion relation is simplified as 


Wo 2k ( Ip 2dk\ 2 {d + lp)lp 4 

At the upper critical velocity 


(25.30) 


K _ 2_ ~ _ 1 _ 

fn,a dk 

the maximum value of oj, vanishes: w, = dui,Jdq = 0. For the velocity higher than 
14, u) q always stays negative. The critical wavenumber is given by 


/ 


1 + 2 k 
dkl T ’ 


(25.31) 


g a oc [fc(l + 2fc)d 3 lx] X! 


(25.32) 
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25.3 Eckhaus Instability 

When the interface deforms sinusoidally as z = ((x,t) + C(x,t) = f Re(a q e' QX )dq, the 
Fourier transform a q (t) develops as 


da q 

dt 


— OJqClq 


(25.33) 


in the linear approximation. Here u q is determined from the dispersion relation 
(25.26). When the pulling velocity V is a little larger than the critical value V c , 
the interface is unstable for the sinusoidal modes with u q > 0. Around the critical 
wavelength q c , u q is expanded as 


U) q = u c - C(q - q r ) 2 


(25.34) 


with u) c = u) qc , and C = —^(d 2 u/dq 2 ) qc is positive. u c is proportional to V - V c . From 
Eq.(25.33) and (25.34), the time evolution of the interface ((x, t) can be approximated 
as 


3£ (x,t) 
dt 


= u r ( — C 



(25.35) 


Near the critical point V c , the unstable mode has the wave numbers around q c , and 
the deformation is expressed as 


<^=A{x,t)^ x (25.36) 

with the slowly varying complex amplitude A(x, t). It satisfies the linear equation 

As Uc is positive, the deformation A increases, but then the nonlinearity should come 
into play. Since the system is invariant by the transversal translation x —> x + <p 
and the space inversion x —* —x, one gets the nonlinear amplitude equation of the 
Landau-Ginzburg type, 

DA B 2 A 

— = u)r A + c —-a l \A\ 2 A. (25.38) 

There is a systematic derivation of the Landau coefficient ai by the reductive per¬ 
turbation method [199]. To specify the meaning of the Landau coefficient aj, we 
consider the amplitude of the critical mode: A(x,t.) — A(t). If ot\ > 0 the third 
order term acts as to limit the amplitude A. The amplitude increases gradually as 
A ~ u'J 2 ~ {V — V c yt 2 from zero near the critical point, as shown in Fig.25.6. The 
transition is similar to the second order phase transition in equilibrium case, and is 
called supercritical. If a\ < 0, one needs still higher order terms to obtain a finite 
amplitude of deformation. For example, with the fifth-order term A 5 with negative 
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Figure 25.6: Schematics of the modulation amplitude of the interface near the insta¬ 
bility. For a positive Landau coefficient an, the bifurcation from the planar to the 
modulated interface is supercritical, whereas for a negative it is subcritical. 


coefficient, the amplitude jumps to a finite value at V c with the hysteresis. The tran¬ 
sition is similar to the first-order phase transition in equilibrium, and the bifurcation 
is called subcritical. We consider hereafter only the supercritical case with a t > 0. 

One can easily find that there is a stationary solution to Eq.(25.38) in the form 
of a single sinusoidal mode as 

A 0 {x) = (25.39) 


Its amplitude A q is determined as 


A 2 _ Ur. - C{q - Qc) 2 = UJq 
q Qi «1 ’ 


(25.40) 


In the domain of w q > 0 where the flat interface is unstable, a stationary profile 
of the interface is possible which is deformed sinusoidally as 2 — Rc(A q c i,lx ) with q 
dependent amplitude A q . But what periodicity is selected? Are all these nonlinear 
solutions stable? 

We consider now the linear stability of the stationary state against the modifica¬ 
tion: 

A{x,t) = [1 + £(*,*)] (25.41) 

with a complex modification For small £, the modification can be written as 
1 + £ ~ e ( ~ e^e' 1 "^. This means that the real part Re£ represents an amplitude 
modification to A q , and the imaginary part Im( represents a phase modulation. By 
inserting Eq.(25.41) to Eq.(25.38), and taking the first order of £, one gets 


dt 


= C 




- ai^(£ + c.c.). 


(25.42) 


This linear equation can be analyzed by using the Fourier transformation as 
£(x,t) — e n< (B x c iQx + B^ iQx ) 


(25.43) 
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with real B\ and Bo. By comparing coefficients for e'^* and c ,Qx , we get the linear 
equation for B\ and B 2 as 

QB l = C [-2 Q(q - q c ) - Q 2 } 5, - [u e - C(q - q c ) 2 ] (B l + B 2 ) 

QB 2 = c[2Q(q-q c )-Q 2 }B 2 -[uJ,-C(q-q l .) 2 ](B l +B 2 ). (25.44) 


Here ai A 2 is replaced by the corresponding term by Eq.(25.40). On requiring that 
Bi and B 2 has nontrivial solution, the eigenvalue equation 


0 -)- 2 CQ(q - q c ) + CQ 2 + uj c - C(q - q c ) 2 u c - C(q - q c ) 2 
u) c - C(q - q,.) 2 U - 2CQ(q - q c ) + CQ 2 + u c - C(q - q c ) 2 


= ft 2 + 2 C[Q 2 + Q 2 m 
= 0 


(q - 9 ,) 2 ] n + 2 C 2 Q 2 


Ql - 3 Q 2 + 



(25.45) 


is derived. Here Ql = oj c /C. Eigenvalues are obtained as 

n± - -C [Ql -(q- qc) 2 + Q 2 ] ± C^/[Q 2 m -(q-q c ) 2 } 2 + 4Q 2 (q- qc y. (25.46) 

For a long wavelength modulation with a small wavenumber Q, one gets 

« -2 C [Ql -(q- g c ) 2 ] + 0{Q 2 ) « -2u q 

O ~ p Q'm ~ 3(? ~ Qc) roc 4 y\ 

" + ~ C Q (26 - 47> 

with corresponding cigenmodes; B 2 _ = B\- = B_ and B 2 ,+ = —Bi, + = B + , or 
= B^(P~ l cosQx and £ + = iB + e n+t sin Qx. Therefore, - mode corresponds to 
the amplitude modulation and + mode corresponds to the phase modulation. Since 
< 0 for u) q > 0, amplitude modulation always damps out. On the contrary, fl + 
can be positive for those modes with wavenumbers q in the region 


Ql>(q-qc) 2 >^f = ^- (25.48) 

The phase modulation destabilizes the stationary solution. Within the region (25.48) 
near the neutral line q = ^(K), the phase diffusion coefficient is negative and the 
phase diffusion mode increases. This is called the Eckhaus instability [55]. Neutral 
curve q = qN{V) is determined from the condition w, = 0. With Eq.(25.34), then, 
<7fv,± = Q< ± \fu) c /C. From Eq.(25.48) the Eckhaus boundary is obtained as qE,± — 
q c ± -j^(qN,± ~ Qc), as depicted in Fig.25.7. 

In the analysis the dispersion relation, Eq.(25.34), is assumed to have a maximum 
at the critical wavenumber q c even in the supercritical regime V > V c . In reality, 
the most unstable wavenumber q m shifts rapidly depending on V. Therefore, the 
approximate Eckhaus boundary can be written around the most unstable mode as 
q E ,± = <i m ± - </,„)• Brattkus and Misbah calculated the Eckhaus boundary 

more precisely, which is quite different from parabola [29]. 
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Figure 25.7: The sinusoidal modulation with wavenumbers q inside the region E is 
unstable against the phase diffusion, and this is called Eckhaus instability. 

25.4 Fully Nonlinear Behavior 

Near the critical point where the surface deformation is small, the amplitude equation 
(25.38) with cti > 0 describes the crystal shape correctly. When the Landau coeffi¬ 
cient ai is negative, or the system is far from the critical point, the deformation is no 
more small and the full nonlinearity has to be considered. In such a case, numerical 
simulation is appropriate to study the large interface deformation [183, 106, 167]. 
We extend the simulation algorithm explained in Appendix A24.3 to the periodically 
modulated surface structures under a constant temperature gradient. For a small ve¬ 
locity near the critical velocity V c , the cellular structure appears as shown in Fig.25.8a. 
On increasing the pulling velocity, the cell groove deepens and the surface takes the 
form of cusp arrays, as in Fig.25.8b. Further increase of the pulling velocity sharpens 
the tip and the dendrite array is formed as in Fig.25.8c. Here the dendrite structure is 
controlled by the anisotropy of the surface tension and the tip radius t and the velocity 
satisfies the solvability condition (24.5). The simulation results qualitatively agrees 
with the experimental sequences of morphology changes as shown in Fig.25.9. 





Figure 25.8: Interface structures in unidirectional solidification obtained by simula¬ 
tion: (a) arrays of cells (V/V c = 1.01), (b) of cusps (V/V c = 1.15), and (c) of dendrites 
(V/V c = 17.6) 







Figure 25.9: Interface structures in unidirectional solidification of snccinonitrile in 
acetone, (a) Arrays of cells, (b) of cusps and (c) of dendrites [180]. 



Figure 25.10: Tilted unidirectional solidification in the pivalic acid-ethanol system. 
Crystal axis istp = 40.5° . The velocity in f.tm/s are (a) 0.5, (b) 1.0, (c) 2.75 and (d) 
10.0 [182], 
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Figure 25.11: Tilted unidirectional solidification at various pulling velocities obtained 
by simulation. The crystal axis is ip = 17° off from the direction of the temperature 
gradient. The tilt angle <p of the tip profile is (a) 4.3°, (b) 8.7°, (c) 12.5°, (d) 14.2° 
and (e) 16.7° [148]. 

In directional solidification there are two tip stabilizing effects, the thermal gradi¬ 
ent and the surface stiffness. What will happen when two anisotropies prefer different 
orientations: the temperature gradient forces the crystal to grow in the z direction, 
but the capillarity prefers another direction with a tilt angle ip. It means that the 
capillary length depends on the orientation ip of crystal axis, and is written as 

d — do [1 - e cos 4(6 - ip)] (25.49) 

at the point with the angle 6 of the interface normal to z axis. 

Tilted unidirectional solidification is observed experimentally [73, 182, 27, 2,144], 
Figure 25.10 shows interface structures at various pulling rates for the pivalic acid- 
ethanol system [182]. On increasing the pulling rate, the interface tilts from the 
direction of the temperature gradient, and approaches to that of the crystalline axis. 
Since the tilting increases by increasing the velocity, the surface kinetics is supposed 
to be relevant. Also from the linear stability analysis, the anisotropy in the surface 
stiffness is found to be irrelevant for the tilting [51, 202], But the recent numerical 
simulations of the unidirectional solidification with local equilibrium assumption [2, 
148] show that the tilting is possible when the crystalline axis is off-angled from the 
pulling direction. The simulation in Fig.25.11 shows that for small pulling velocity 
V the cell tip is oriented close to z axis, the orientation of temperature gradient, 
but for large V the dendritic crystal is oriented close to the crystalline axis. Thus 
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the experimental tendency of large tilt angle at fast pulling can be reproduced by 
the surface stiffness, too. More studies are necessary to identify the real mechanism 
of crystal tilting, whether it is due to the anisotropy in capillarity or to that in the 
kinetics. 


26 Eutectic Growth 

For some alloy system such as Pb-Sn, the phase diagram looks as shown in Fig.26.1. 
Due to the mixing entropy the melting temperature of an alloy often decreases from 
the pure material. When two liquidus lines of AB alloy decrease down from the pure 
A and from the pure B melting temperatures by mixing the other component, they 
meet at some concentration of B species c E and the temperature T E . At Tr, the 
liquid solution with the concentration c E coexists with two crystals, phase a with 
concentration c|(< c E ) and phase 3 with concentration ef(> cr). This triple point is 
called the eutectic point. At this eutectic point the melting temperature is minimum. 

If one grows this eutectic alloy unidirectionally in a Hele-Shaw cell, what kind 
of structure appears? Starting from the liquid with the eutectic concentration c E , 
crystallization of a or f) phase alone cannot satisfy the material conservation, a and 
3 phases should appear simultaneously. It may intuitively expected that a and ft 
phases appear alternatively as lamellae parallel to the growth direction: When the 
crystal is growing in positive z directions, the symmetry in x and —x directions yields 
that the a/3 phase boundary is expected to align in 2 direction and perpendicular to 
the liquid-crystal interface in x direction. This is the structure found in experiment 
[91] at low pulling velocity, as shown in Fig.26.2 [108], The problem to be posed is 
the selection of the periodicity A. There is a detailed calculation by Jackson and 


T 



Figure 26.1: Equilibrium phase diagram of an eutectic AB alloy. 
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Figure 26.2: Lamellar growth of eutectic crystal, CBr 4 - C 2 C1 6 [108] 


Hunt, which is reproduced in Appendix A26 [91]. Here a phenomenological discussion 
is given on the periodicity selection. 


26.1 Lamellar Structure 

Since the thermal diffusivity is larger than the chemical diffusivity, the temperature 
gradient Gt = dT/dz can be assumed constant. By choosing the origin of z axis at the 
point where the temperature is at the eutectic value T E , the temperature distribution 
is described as 

T{z) = T e + G t z. (26.1) 

The liquid (z > 0) is hot and the crystal (z < 0) is cold. We consider the simplest 
case that a crystal is growing from a solution with eutectic concentration c E . When 
the crystal takes a lamellar structure with a periodicity A, the ratio of a phase r] and 
that of 0 phase 1 — 77 should satisfy the conservation relation 

ce = f?c<- + (1 - rj)4 ( 26 - 2 ) 


or 


with the miscibility gap 



Ac = cf - eg. 


(26.3) 

(26.4) 


Since eg is smaller than c E , B atoms are expelled out from the a crystal into the liquid 
by crystallization, and the liquid concentration c“ in front of a phase increases above 
c E . On the contrary, in front of the 0 phase, B atoms are sucked up by 0 crystal and 
the liquid concentration decreases below c E during crystallization. As is apparent 
from the magnified phase diagram Fig.26.3 with metastable liquid branches, deviation 
of the liquid concentration from the eutectic value means that the supercooling at the 
interface under the local equilibrium assumption. The crystal concentration eg and 
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Figure 26.3: Phase diagram near the eutectic temperature T F: with metastable 
branches of coexistence curve. 

correspondingly differ from those at the eutectic temperature, but their differences 
are small and do not modify the conclusion essentially. 

Since the concentration in front of each crystal is different from that at z —> oo 
or Ce, the concentration diffusion takes place. The diffusion flow compensates the 
material deficit or surplus produced by the crystallization. In front of the a crystal, 
v n {dl - cjj) « v n (c F - eg) of B atoms are expelled from the unit surface per unit time. 
This excess is transported by the diffusion flow -D c dc/dn. However, the material 
does not need to be transported up to z -> oo: It is sufficient to be transported 
only to the neighboring lamellae of (i phase, where the material is deficient. Since 
the concentration difference in liquid in front of the a and (j phases is eg — of for a 
separation of A/2, the growth rate is determined as 

(c E -c%)v^Q a D r ^f. (26.5) 

If the interface is flat, the direction of material diffusion x is orthogonal to the growth 
direction z and the diffusion does not contribute to the growth. Actually, the interface 
is curved, and the contribution to growth appears. This is taken into consideration 
in Eq.(26.5) by the factor Q„. Its detailed form is explained in Appendix A26. In 
consideration of the material conservation in front of the f) phase, we get the similar 
formula with coefficient Qp. In the steady state where the growth velocity of a and 
fj phases are the same, the coefficients should satisfy the relation (1 — r))/Q a =s r]/Qp 
due to the conservation (26.2). 

The liquid-crystal interface is curved in order to satisfy the mechanical force bal¬ 
ance at the triple point, where three phases, liquid, a crystal, and ft crystal phases, 
meet (Fig.26.4). The force balance in x and z directions is described as: 

7 i.r» sin 9 a + 7 L/j sin 9 P = i af) (26.6a) 

7i, a cos 6 a = 7r,/? cos 9g. (26.6b) 

Here 7 Lq is the surface tension between the liquid and the a phases which is assumed 
isotropic for simplicity, and so on. In equilibrium the liquid- a (ft) crystal interface is 
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Figure 26.4: Interface profile near the triple point for the eutectic growth. 


part of a circle with a radius R a (Rg), and the angle 0 n ( f)g ) is related to R n (Rg ) 
and the periodicity A as 

2R q sin 0 a = r)\, 2Rpsin0p = (1 — r/)A, (26.7) 


as is evident from Fig.26.4. 

The interface temperature is given by the concentration deviation from the eutec¬ 
tic values and also by the Gibbs-Thomson curvature effect as 

T a = Te — m a (C], - ce ) - (26.8a) 

U a 

^ = r E + m^- CE )-^-i-. (26.8b) 

Lp Up 

Here m a = |dT a /dc£| and mp = | dTp/dc^\ are the absolute slopes of the liquidus lines 
for a and (5 phases, L a , Lp are latent heats of a and /3 phases, respectively. 

By adding Eq.(26.8a) and (26.8b), liquid concentration difference is written as 


cf -cf = + 

\m a m p J 


/ IX.gTp, 1 ^ t^Tr 1 \ 

\m a L a R a mpLpRp)' 


(26.9) 


Here the interface supercoolings are defined as AT" = T E - T" and AT 1 * = Tr — T fi . 

We consider the simplest case where crystal phases a and ft are symmetric: Cf, = 
1/2, m Q = mp = m, 6 a = Op - 0, r) = 1/2, Q„ = Qp = Q , 7l„ = 7i./i = 7-. AT“ = 
AT' 9 = AT, and c£ - eg = Ac/2. By using the capillary length d = 'yx.JTyJLmAc., 
the interface supercooling is written from Eqs.(26.5) and (26.9) as 


AT = Tp, 



(26.10) 


which consists of two contributions; Diffusional 011 c ATd = T Vi a n X(l n and the kinetic 
one AT k = T^and/X. Here In = 2 D c /v is the diffusion length, and material parame¬ 
ters are represented as ap = mAc/4QT F ,, Ok = 4 sin 0mAc/T K . The relation between 
AT and A is shown in Fig.26.5. 
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Figure 26.5: Minimum undercooling hypothesis for the selection of the eutectic lamel¬ 
lar periodicity. 

In order to determine the lamellar periodicity A, the minimum supercooling hy¬ 
pothesis is used by Jackson and Hunt: 3A TjdX — 0. Then the periodicity is deter¬ 
mined as 


This yields the relation between the pulling velocity v and the periodicity A m as 

vX - = wte =consL (2612) 

Also the interface undercooling A T m is determined as 


(26.11) 



Figure 26.6: Variation of the average interlamcllar spacing A with the growth velocity 
v for the carbon tetrabromide (CBr 4 )- hexachlorocthenc (C 2 C1 6 ) eutectic system [174], 
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A T m = 2T EV 'a D a K 


d_ 

Id 


(26.13) 


and it satisfies the relation with the pulling velocity v as 


AT* 


v 


2 sin# 

~~Q~ 


(mAc) 2 — = const. 


(26.14) 


The full nonlinear analysis by computer simulation [107, 95, 98, 96] also show 
that the interface undercooling has the minimum as a function of the periodicity, and 
that the Jackson and Hunt theory provides the accurate position of the minimum 
undercooling [96]. The scaling relations, (26.12) and (26.14), axe also shown there 
[98, 96]. These relations are also confirmed by experiments, as shown in Fig.26.6 
[174, 203]. 


26.2 Parity Breaking and Oscillation 

On increasing the pulling velocity V or on varying the liquid concentration from the 
eutectic value, the lamellar structure is found to be modified. 

For slow pulling rate, the phase boundary between a and /? phases align parallel 
to the temperature gradient. On suddenly increasing the pulling rate by a factor 4, 
the periodicity selected by the minimum undercooling condition should be half of the 
initial one, as Eq.(26.12) tells. But such a large structural variation is not possible 
topologically. The numerical simulation using the boundary element method in the 
stationary code found that the parity (left-right symmetry) will break for large pulling 
rate [52,95, 97]: The a(3 or 0a triple point shifts transversally along the liquid-crystal 
interface, and the a(3 phase boundary tilt from the temperature gradient, as shown 
in Fig.26.7a. The tilting is found later in many experiments (Fig.26.7b) [57, 58]. 



(a) (b) 


Figure 26.7: Parity broken lamellae in eutectics obtained (a) by simulation and [97] 
(b) by experiment of CBr 4 —C 2 CI 6 [57]. 
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Figure 26.8: Periodic oscillation of lamellae in off-eutectic alloy (a) by the experiment 
of AI-CUAI 2 , and (b) by simulation [204], 

This tilting is different from that discussed in the previous section 25.3. There the 
tilting is enforced by the tilting of crystalline axis from the growth direction and the 
parity is externally broken. Here the system has the left-right parity symmetry, but 
the realized state breaks the parity spontaneously. 

By shifting the. far field concentration c M from the eutectic value c®, system shows 
the nonstationary structure, as periodic oscillation shown in Fig.26.8a [204], The 
linear stability analysis showed the instability of the oscillatory motion of the triple 
point [53, 94]. There is a model simulation which shows the oscillatory behavior, and 
the result is in good agreement to the experiment as shown in Fig.26.8b [204], 


27 Diffusion Effect on Polyhedral Crystal: Berg 
Effect 

So far it is assumed that the surface kinetics is infinitely fast and the local equilibrium 
is realized at the interface in order to stress the diffusion effect on the front instability. 
This approximation is not far from reality for an ice crystal growing in water [65, 66]. 
The prism face of ice in water is rough, and the heat transport controls the growth. 
On the other hand, for snow growing from the water vapor, the approximation does 
not hold since the surface of the snow is sharp and faceted. The kinetics is slow and 
has to be properly taken into account. 

If the interface is atomically rough, the kinetic coefficient K in Eq.(18.8) is finite 
but not singular. The growth shape can be anisotropic but is nonsingular and smooth. 
The orientation dependence of K leads to the new scaling relation of the dendritic 
growth rate and the tip radius [32, 33, 50, 169], When the preferred orientations of 
the surface tension and of the kinetic coefficient are different, rich variety of growth 
shapes with dynamical shape transitions is expected [45] 
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Figure 27.1: (a)Electrostatic potential around metallic polygon. The metal surface is 
equipotential. (b) Concentration field distribution around the square crystal growing 
steadily with a constant normal gradient q n = dc/dn along the surface. 

If the interface is atomically flat, the kinetic coefficient K is singular, and the 
crystal is expected to take polyhedral morphology covered with flat facets of singular 
surfaces. The growth of these singular surfaces is governed by the surface kinetics 
discussed in Part III. However, material has to be transported to the growing interface 
by the diffusion, and the diffusion also influences the growth form. One has to consider 
both the kinetic and diffusion processes in this situation. 

First we consider the concentration distribution around the polyhedral crystal 
growing from the solution. It is determined by Eq.(18.6-18.8). Usually, crystal growth 
rate is small for kinetic-controlled growth, and one can assume the stationary condi¬ 
tion that the diffusion field quickly adjust its distribution around the given crystal 
morphology: du/dt — 0. Thus the diffusion field u satisfies the Laplace equation 

V 2 m = 0 (27.1) 

instead of the diffusion equation (18.6). An analogy holds between the concentra¬ 
tion distribution «(r) and the electrostatic potential 0(r). Far from the crystal, 
the concentration distribution reflects the isotropic nature of the system and equi- 
concentrations are asymptotically spherical. If the concentration is assumed constant 
on the polyhedral face, u(r) corresponds to electrostatic potential <f>( r) around the 
metallic polyhedron, and looks as shown in Fig.27.1a. Near the corner of the polyhe¬ 
dron, the spacing between consecutive equi-concentrations are narrow, and the normal 
growth rate determined by Eq.(18.7) or V n = -D c dc/dn is larger at the corner than 
at the center. Then the polyhedral face cannot remain flat. To keep the polyhedral 
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Figure 27.2: Concentration distribution around the growing NaCl crystal from water 
solution [24]. 

face flat, there should be concentration difference on the surface, especially between 
those at the corner and at the center. 

If we assume the constant normal gradient on a flat but bounded surface, the 
concentration at the corner should be higher than that at the center, as shown in 
Fig.27.lb. This was shown analytically in two dimensions [173], Experimentally, the 
concentration profile around the polyhedral crystal is observed by Berg as shown in 
Fig.27.2, and the concentration difference along the surface is observed. This effect 
is called the Berg effect [24]. The distribution of the supersaturation along the flat 
surface is also measured [42]. 

If the face is really flat and the concentration varies on it, the kinetically controlled 
growth rate (18.8) varies from position to position, and it is impossible to keep the face 
flat. To realize a steady growth, the kinetic coefficient should vary along the surface. 
This is possible when a macroscopically flat surface contains many microscopic steps 
and is in fact a vicinal face. Then the kinetic coefficient depends on the concentration. 
For example, at the corner of the crystal where the supersaturation is high, the two- 



Figure 27.3: Formation of snow dendrite by simulation [201], 
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dimensional nucleation supplies steps toward the center of the face [119]. The step 
advancement velocity v decreases near the center where the supersaturation is low, 
but if the step density is high and the step separation l is small, the normal growth 
rate R = avjl can be kept constant over the surface. 

If the supersaturation at the corner is too high or the crystal has grown too 
large, the difference in supersaturations at the corner and at the center becomes too 
large to keep the constant normal growth rate over the flat face. The corner starts 
to grow faster than the center, and the instability sets in [47, 119]. This initiation 
of the corner instability is analyzed numerically, and the corner instability and its 
size dependence is found, if the nucleation controls the kinetics [119]. Also by using 
the boundary element method the growth of a kinetically-controlled two-dimensional 
crystal is simulated. The circular crystal becomes polygonal, and then produces the 
faceted arm from the corner of polygon, similar to the snow growth, as shown in 
Fig.27.3 [201]. For the sidebranch formation of a faceted dendrite, however, more 
studies seem necessary. 



Part V 

Appendices 

A9 More on Surface Roughening 
A9.1 Solid-on-solid (SOS) model 

In order to focus on the temperature dependence of the surface structure, the solid-on- 
solid (SOS) model is often used [112, 43], The model picks up the freedom of surface 
height z(i) such that no vacancy in the crystal, no overhangs at the surface and no 
crystal clusters in the vapor are allowed. Here i is the site in a d-dimensional hyper- 
cubic lattice with 2 d nearest neighbors, and the crystal occupies d + 1 dimensional 
space. In a unit of an atomic size a, the height h(i) = z(i)/a takes integer values, 
h(i) = —oo, • • •, — 1,0,1 • • •, oo. Since the height difference between neighboring sites 
costs energy, the Hamiltonian H is written as 

H = J'Z\Ki)-h(j)\”. (A9.1) 

<ij> 

Here the model with p = 1 is called the (absolute) SOS model, the model with p = 2 
the discrete Gaussian (DG) model. The latter is convenient to treat analytically. If 
the height takes only two values, h(i) = 0,+1, the model corresponds to the Ising 
model. 

As to measure the surface roughness, the height difference correlation function 
[49] 

G(ry) = i(h(i) - h(j)) 2 ) (A9.2) 

is used. When the surface is flat, heights arc strongly correlated and G remains finite 
even for r,j —> oo. When the surface is rough, heights at different positions fluctuate 
independently, and G(r —> oo) diverges. In an Ising model, G can at most be 1, and 
the surface can never be rough in this sense. 

A9.2 Monte Carlo simulation of the SOS model 

This appendix contains a source list of an example program of Monte Carlo simulation 
explained in Section 9.1. It is hoped that it provides a first step for the reader to 
get familiar with the simulation and to develop one’s own program. Random number 
generator is copied from the reference [67]. 
program sos 

c Monte Carlo simulation of 2D-sos model for roughening 
c============================================= 

parameter (lxmx=128,lymx=128) 
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dimension ih(lxmx,lymx),prob(6) 
read(*,*) lx,ly,loop,temp,field,iseed 

c- 

c ih(ix,iy); solid height 

c lxmx.lymx; maximum linear size in x- and y-direction 
c lx.ly ; actual linear size of the system 
c prob(k); transition probability 
c k=i for crystallization 

c k=2-6 for desorption with k-2 neighbours 

c temp; temperature, (coupling J=l) 
c field; chemical potential difference 

c loop; loop number for small accumulation and averaging 
c iseed; seed of the random number 
c mag; magnitude of height sum 
c nnex; nearest neighbor energy 

c- 

call ranint(iseed) 

call init(lxmx,lymx,lx,ly,ih) 

write(*,*) ’lx,ly,loop,temp,field;’ 

write(*,*) lx,ly,loop,temp,field 

write(*,*) ’ilp,mag,nnex’ 

call engy(lxmx,lymx,lx,ly,ih,mag,nnex) 

call prbini(prob,temp,field) 

c- monte carlo loop 

do 10 ilp=l,loop 
do 11 inlp=l,lx*ly 

call me(lxmx,lymx,lx,ly,ih,prob,mag,nnex) 

11 continue 

write(*,*) v ilp,mag,nnex 
10 continue 

c- check the final height and exchange energy 

call engy(lxmx,lymx,lx,ly,ih,magi,nnexl) 

write(*,*) mag,magi,nnex,nnexl 

stop 

end 

subroutine init(lxmx,lymx,lx,ly,ih) 

c--- 

c initiallization 

c- 

dimension ih(lxmx,lymx) 
do 200 iy=l,ly 
do 200 ix=l,lx 
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ih(ix,iy)=0 
200 continue 
return 
end 

c=================================—================ 

subroutine me(lxmx,lymx,lx,ly,ih,prob,mag,nnex) 

c- 

c Monte Carlo step 

c- 

dimension ih(lxmx,lymx),prob(6) 
isum=0 

ix=randm()*lx+l 

iy=randm()*ly+l 

ixp=ix+l 

ixm=ix-l 

if(ixp.gt.lx) ixp=ixp-lx 
if(ixm.lt.1) ixm=ixm+lx 
iyp=iy+l 
iym=iy-l 

if(iyp.gt.ly) iyp=iyp-ly 
if(iym.lt.l) iym=iym+ly 
if(ih(ixp,iy).It.ih(ix,iy)) isum=isum+l 
if(ih(ix.iyp).It.ih(ix,iy)) isum=isum+l 
if(ih(ixm,iy).lt.ih(ix,iy)) isum=isum+l 
ifCihCix,iym).It.ih(ix,iy)) isum=isum+l 
trpr=prob(1)+prob(isum+2) 
trtr=randm() 

c- adsorption 

if(trtr.lt.probCl)) then 

ih(ix,iy)=ih(ix,iy)+l 

jbup=0 

if(ih(ixp,iy).It.ih(ix,iy)) jbup=jbup+l 
if(ih(ix,iyp),lt.ih(ix,iy)) jbup=jbup+l 
if(ih(ixm,iy).It.ih(ix,iy)) jbup=jbup+l 
if(ih(ix,iym).It.ih(ix,iy)) jbup=jbup+l 
mag=mag+l 

nnex=nnex+2*j bup-4 

return 

endif 

c- desorption 

if(trtr.lt.trpr) then 
ih(ix,iy)=ih(ix,iy)-l 
mag=mag-l 
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nnex=nnex+4-2*isum 

return 

endif 

c - otherwise 

return 

end 


subroutine engy(lxmx,lymx,lx,ly,ih,mag,nnex) 


c energy calculation 


dimension ih(lxmx,lymx) 

mag=0 

nnex=0 

do 100 ix=l,lx 
ixp=ix+l 

if(ixp.gt.lx) ixp=ixp-lx 
do 100 iy=l,ly 
iyp=iy+l 

if(iyp.gt.ly) iyp=iyp-ly 
mag=mag+ih(ix,iy) 

nnex=nnex+iabs(ih(ix,iy)-ih(ix.iyp))+iabs(ih(ix,iy)“ih(ixp,iy)) 
100 continue 
return 
end 


subroutine prbini(prob,temp,field) 
dimension prob(6) 
ef=exp(field/temp)+exp(4,/temp) 
prob(l)=exp(field/temp)/ef 
do 10 nn=2,6 

prob(nn)=exp(2.*(nn-4.)/temp)/ef 
10 continue 
return 
end 

subroutine ranint(ix) 
common /rand/ m,j 
dimension m(521),ia(521) 
do 10 i=1,521 
ix=69069*ix 
10 ia(i)=sign(l,ix) 
do 20 j=l,521 
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ip=mod((j-l)*32,521)+l 
m(j)=0 
do 30 i=l,31 
ii=mod(ip+i-2,521)+1 
m(j)=2*m(j)+(ia(ii)-l)/(-2) 
ij=mod(ii+488,521)+l 
30 ia(ii)=ia(ii)*ia(ij) 
ii=mod(ip+30,521)+l 
ij=mod(ii+488,S21)+l 
20 ia(ii)=ia(ii)*ia(ij) 
j=0 

return 

end 

function randmO 
common /rand/ m,j 
dimension m(521) 
ip=521 
iq=489 

j=j + l 

if(j.gt.ip) j=l 
k=j+iq 

if(k.gt.ip) k=k-ip 
m(j)=ieor(m(k),m(j)) 
randm=m(j)*0.4656612e-9 
return 
end 


A9.3 Continuous Gaussian model 

If the height variable h(i) does not take the discrete integer values but takes contin¬ 
uous values, and the energy is given by (A9.1) with p = 2, the model is called the 
continuous Gaussian model. The partition function and correlation function of the 
model can be calculated straightforwardly. 

The partition function is defined as 


Z = n'I 1 r dh(i) exp 

•J —oo 


- J £ 

<o> 


(h(i) - h{j)) 2 


k R T 


One can introduce the Fourier transformation as 

1 


h(i) = y^Ee - ’ ,r, 7i(9) 


(A9.3) 


(A9.4) 
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with wavenumbers q = 2tt n/L and n = -(L - l)/2, • • • ,0,1, • • ■, L/ 2. Here L is the 
linear dimension of the system such that N = L d . The coefficient h(q) is given by 

h(q) = ~^ <qri Ki). (A9.5) 


In terms of h(q), the partition function is written as 


z = n, jT dft(fl) exp[-i£ = nj2nGo(q) 


(A9.6) 


Here C?o(<j) is the lattice Green’s function defined by 


Go (q) = 


k B T 1 __ k B T 1 

2J E«(l - «**«) * 2J 


(A9.7) 


where $ is a vector connecting nearest neighbor sites. Surface free energy is obtained 
as 


1 Nr ri d n 

F = -2*bT21ii[2»Go(9)] = ~jk B T j —^ln27r G 0 (q). (A9.8) 


(2?r) d 


In the last equality, the summation is replaced by the integration for a large sys¬ 
tem size N as = N d d q/(2n) d . The height difference correlation function is 
calculated for large r as [115] 


G(r) = ([h(r) - h(0)) 2 ) = |r£(|%)| 2 )(l -cos?r) =2 J^G 0 (g)(l 


k B T 


1 


r 2-d 


1 


J 2 d_ 1 7 r d / 2 r(d/ 2 ) 2 — d 
k B T 

-—- In r + const 
2irJ 


+ const + 0(r‘ d ) 


for d / 2 
for d — 2 , 


cos qr) 
(A9.9) 


where T(d/2) is the Gamma function [1], The d-dimensional surface with d < 2 
is always rough because the height difference correlation function diverges for two 
separate points, r —* oo. On the contrary, for d > 2 the surface is always smooth. 


A9.4 Variation approximation 

Continuous Gauss model shows that the surface is rough in two dimensions. In the 
real crystal with d = 2, the surface can be both smooth and rough. The discrepancy 
is due to the neglect of the discreteness in the surface position in the continuous 
Gaussian model. The height should be quantized in unit of atomic size a. For 
analytical treatment on the phase transition and critical phenomena, a modification 
is introduced in the model Hamiltonian as 

n = jJ2 [*(*) - - i/E C0S ( 27rft (*)), 

<ij> i 


(A9.10) 
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Figure A9.1: Convexity of exponential function. 


where the height variables h(i) are continuous, but the additional one-body potential 
with a small positive coefficient U prefers the height to take integer values. 

We here study the phase transition of the model (A9.10) by means of the varia¬ 
tion method [59]. The theory depends on the convexity of the exponential function; 
e -(a+6 )/2 < (e _ “-j-e _ft )/2. This inequality is obvious from Fig.A9.1. When the stochas¬ 
tic variable x is distributed according to the probability Pr(or), the average (x) is given 
by ( x) = JPr(x)xdx, and the inequality is generalized to < (e~ x ). 

For a thermodynamic system with a Hamiltonian Ho, the canonical probability is 
defined by 

Pr = Zq 1 exp ( — (A9.ll) 

with 

z °=' Ii “ p (-j£f) - exp (-s?) - (A9 - 12) 

When the true system is described by the Hamiltonian 

H = H 0 + V, (A9.13) 


one can show that 


exp 



s <» p (-^)>—(£h» p (- 

= exp (tS : )“ p (^ts i )^ 


Ho + 
k B T ) 

(A9.14) 


Therefore, the true free energy F satisfies the following inequality 


F<F 0 + (V) 0 = F*. 


(A9.15) 


If we make an appropriate choice of Ho which minimizes F*, then it should be a 
good approximation to the true free energy F. In order to make the configuration 
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summation tractable, we choose the effective Hamiltonian Ho as 


Wo = \k B T^ 

Z Q 


mi 

G(q) ' 


(A9.16) 


where q —dependent parameters G(q) are to be determined to find a minimum F*. 
With this choice of Ho, the probability Pr given by Eq.(A9.11) is Gaussian with 
(h(q)) = 0 and <|/i(g)| 2 ) = G(q). 

Fq is calculated straightforwardly and is equal to (A9.8) with the replacement of 
Gq( q) by G(q). (Ho)o is simply Nk B T/ 2. In (H) o, the average of the first term is 
easily calculated as 


k B T 


£ um - Mi)) 2 ) = ~££<IM?)| 2 )(1 - e"*) = 1^. (A9.17) 


<t}> 


Go{q) 


Since there are only first and second order cumulants for a Gaussian distribution, the 
term containing the one-body potential can be calculated by using the relation 


(cos2frh) 0 = Re(e 2,r,h )o = exp 


i2ir{h)o - ~-(h 2 )o 


= exp [-27T 2 (h 2 ) 0 ] (A9.18) 


with 


(M)o = ^£(IMI 2 )o = ^£G( 9 ). 

q q 


The approximate free energy F* is obtained as 


(A9.19) 


F* 

k B T 


F 0 + {{H - Ho))o 
k e T 

1 v-, r„ /~*/ 1 v- G{q) NU 

-2? ln l 2 ' G ( ,)l + 5?5SS-^ 


exp 



N 
2 ' 

(A9.20) 


If one calculates the higher order term {V n ) 0 , it is proportional to U n , and the present 
approximation is good for a small perturbation U. To minimize, F* is differentiated 
by G{q) as 


d{F‘/k B T) 

dG(q) 


1 1 


2 n 2 U 

2G(q)^2G 0 (q) + k B T 


1 

• + x 


exp 


O IT 2 

* i 


= 0. (A9.21) 


Since the last term is independent of the wave number q, the ^-dependent coupling 
G(q) has the form 


G(q) = 


K 


Goiq)-' + K~'t 2 q 2 + f~ 


(A9.22) 
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Figure A9.2: Determination of the correlation length £ or the parameter x = (n/£) 2 
at various temperatures t = nknT/AJ. 


where K 1 = 2 J/k B T, and the correlation length £ is determined as 


K~ l C 2 


(2n) 2 U 

k B T 


exp 



ct 2 q K 
(27r) 2 q 1 + 



(A9.23) 


By denoting x = (i r£) 2 , t = nK/2 = -Kk B T/AJ, Eq.(A9.23) is transformed as 


x= 2 -j{Th) t " fix] - (A9 - 24) 

Two curves y = x and y = f(x) are depicted in Fig.A9.2 at various temperatures t. 
The intersection of these two curves gives the solution of Eq.(A9.24). The slope of 
f{x) at x = 0 changes drastically at t = 1, since 


df(x) = 2U t (1 _ 1 \ /_£_y 

dx J \a; 1 + xJ \1 + x) 


x 


i-i 


(A9.25) 


For small U and t > 1 there is only a trivial solution x — 0, but for t < 1 there is 
another nontrivial solution near x ~ 0. For small x, the solution is approximated as 


x = 




V2 U/ 


exp 


\n(J/2U) 
1 -t 


(A9.26) 


For t < 1 the free energy F" for x / 0 is smaller than that for x — 0, and the state 
with a finite £ is stabler than the state with £ = oo. Therefore, the phase transition 
takes place at t = 1 or 

k B T R = — J. 

7T 


(A9.27) 
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The height difference correlation function is calculated by 

G(r) = 2 f -~G{q){\ - e^) « - f A 

J (2?r) 7T ,/miii(jr/fyr/£) ^ ^ 

/"M 


•ln[min(r,0] • 


(A9.28) 


Above the roughening temperature Tr = AJ/irka or t > 1 , there is only a solution 
x = 0 ( or £ = oo), and the height fluctuation diverges at r —* oo: The surface is 
rough. On the other hand, for T < T R or t < 1, f is finite and the height correlation 
function saturates as 


G(r »£) ~ |~j ln £ ~ (1 _< ) 


(A9.29) 


for r —* oo. 

Below Tr the surface height fluctuates within a range of £. This means that the 
steps of a length of order £ can be thermally generated. The step energy (5 is then of 
the order k B T/£ as 

P ~ ~ exp ~ T ^ Z _ T • (A9.30) 

It approaches zero as well as its all derivatives by T on approaching the roughening 
point T —» Tr: It has an essential singularity at T R . For T > Tr the correlation 
length £ = oo and 0 = 0 . 


A9.5 Renormalization group theory on roughening 

In the variation method (Appendix A9.4), only the lowest order of perturbation U 
is considered. The renormalization group method takes the effect of higher order 
perturbation systematically [113, 114, 145, 146]. 

By coarse graining the two-dimensional space, one gets the Hamiltonian in a 
continuum space limit as 

7i = j dxdy ~(Vr) 2 - Vcos^^j = J dxdy -~-(Vh) 2 — Vcos{2irh) . 

(A9.31) 

Here the lattice parameter o is notified explicitly, and the height z(x, y) is normalized 
as h(x, y) = z(x, y)/d. L is the linear dimension of the system. Coupling constants in 
(A9.31) are related to the microscopic parameters in (A9.10) as follows: the surface 
stiffness is here denoted by 7 = 2 J/a 2 and the strength of the one-body potential by 
V = U/a 2 . 

The height h( r) is expanded in Fourier series with wavevectors q between the 
lower, q c = 7 r/£,, and the upper, A = 7 r/a, cutoffs. 

A(r) = £ h,(q)e ,qT = £ e tqr h( q) + £ e ipr h( p) 

q t <q<A q c <q<b^ A 6 -1 A<p<A 

s h<(t) + h> (r). 


(A9.32) 
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In Eq.(A9.32) the height variable is decomposed in the short wavelength component 
h > (r) and the long wavelength component h < (r). The full Hamiltonian consists 
of a term which contains only long wavelength components, that with only short 
wavelength components and that with mixed components as 

H(h{q),h(p)) = H% +?Q +H m = ~ £ q 2 \h{q)f + £ p 2 |/i(p)P 

qc<q<b~ l A ^ b~ l A<p<A 

— V fd,dyco S [2,(h< + h>)]. (A9.33) 

The potential term with the coefficient V is denoted here H m for short. In the 
renormalization method, the short wavelength component h > (r) is integrated out 
and its effect on the long wavelength component h < (r) is systematically studied. The 
renormalized Hamiltonian H < contains only the long wavelength variables h K (r), and 
is defined by 


exp 



= Tr>e- K/A: » r = n, 


6 _1 A<p<A 


/ oo 

dh(p) exp 

-oo 


H(h(q),h{p)) 


k n T 


(A9.34) 


Here Tr > is the trace over the short wavelength component h(p). Actually, the renor¬ 
malized Hamiltonian H < is calculated up to the second order of V. The thermal 
average <• • •} is defined as that by the short wavelength mode with unperturbed 
Hamiltonian Hq as 


Tr> ■ ■ ■ cxp(-Hp/kpT) 
{ "' } Tr > exp( -HZ/k B T) 


(A9.35) 


It is actually a Gaussian average with { h(p )) = 0 and the height correlation (|h(p)| 2 ) = 

k„Thay. 

Average of the term containing the mixed part H m is expanded in the cumulant 
as 




(A9.36) 


The first order of V can be calculated as 


H (1 > 

k B T 


V 

k*T 

V 

k a T 

V 

k B T 

V 


J d 2 r{ cos 27r(/i < + h > )) 

J d 2 r [cos(27rh < )(cos(27rh > )) - sin(27rh < ){sin(27rh' > ))] 

J d 2 rcos(27r/i < )exp [—27r 2 ((h > ) 2 )] 

J d 2 r cos27r/i < (r). (A9.37) 
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Here the effective potential strength V is 

V = V exp [—2 tt 2 ^(0)] = vb- rkBThai (A9.38) 

with the height correlation function of the short wavelength defined by 

~ S^' tl(A,)lrl h - ( A9 - 39 ) 

which decays very fast for r > n/A due to the asymptotic behavior of the Bessel 
function Jo- The second order term ?t (2) is calculated as follows: 

7^(2) l f V \ 2 r 0 /■ „ 

~fc^T = 2 \k^f) J ** T J ^ r> [( cos2nh ‘ c °s2?r/i,> - ( cos2nh}(cos2nhi )] 

« 2 (fT’) I ^(sin27r/i < sm2irhf)(2nh > ■ 2nhf ) 

= J dr J dr { g(r - n) [cos2tt(/i < + hf) -cos2n{h K - hf)], 

(A9.40) 

where the notation h = h(r ) and hi = h(r\) is used. Since g{r - r{) is almost zero 
for |r - n| > 7 t/A, the integrand gives contribution only for ri ~ r. The first term 
produces the higher harmonics cos 2;r(h < + hf) « cos 47 rh < , and is irrelevant for 
further discussion. The second term is expanded as 

J <Pr l Jo(A\r - r l \)cos2n(h < - hf) « const - 4n 3 A ” 4 A j (V/i < ) 2 , (A9.41) 

where A is a complicated but a smooth function in the region of interest [146]. Eqs. 
(A9.40) and (A9.41) shows that the second order perturbation gives the renormaliza¬ 
tion of the surface tension. The final form of the renormalized Hamiltonian has the 
same form with the original one as 

n K = J L dxdy ^(Vh < ) 2 -V r cos(2wh<) 

fL/b 'yn^ _ _ 

= J dxdy — (Vh K ) 2 - V6 2 cos(27r/i < ) , (A9.42) 

where x = bx and y = by. The coupling constant is renormalized from the original 
ones, 7 a 2 and U = V(v/A) 2 , to 

_ 2 2 2A(2/X)U 2 , 

7 a = 7 a 2 + — — In b 

7 a J 


(A9.43) 
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Figure A9.3: The renormalization group flow diagram in X = 2"fa 2 /i:kBT and 
y = 4 U/irksT phase space. The temperature variation of a real system takes place 
along the dashed line. At low temperatures to the right of the separatrix, the po¬ 
tential strength Y increases by renormalization and the surface is smooth. At high 
temperatures to the left of the separatrix, Y vanishes by renormalization and the 
surface is rough. 


U = V = U&-' k » Thal « U + U ^2 - ~Pj In b (A9.44) 

with X — 2')a 2 /nk[ ) T. Therefore, by coarsening the system with the scale e = In 6, 
variables X and Y — W/irkaT arc renormalized as 


dX 

de 

dY 

de 


A{2/X)Y 2 
2 X 



(A9.45a) 

(A9.45b) 


Eqs.( A9.45a,b) are called the renormalization group equations. There is a special fixed 
point X m = 1, y* = 0, and it will be shown that the critical behavior is controlled 
around this point. Therefore, one can assume that the function A is constant with 
its argument at X" = 1 as A(2) = A = 0.398. The flow by renormalization in XY 


plane is obtained by 


dY _ 4 X - 1 
dX ~ A Y 


(A9.46) 


or after integration the renormalization trajectory is obtained by hyperbola as 


l (X -i) 2 -y 2 = -c. 


(A9.47) 


The flow by renormalization is shown in Fig.A9.3. The line Y = 0 for X < 1 is 
a line of fixed points. The separatrix C = 0 or \fAYj2 = 1 — X corresponds a 
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critical curve, and terminates at the end of fixed points at X* = 1, F* = 0. Below 
the separatrix, the strength of the periodic potential U is renormalized to zero, and 
the surface is rough in a large scale. Since the physical system has fixed values of 
surface tension 7 oa 2 and the potential Uo, the temperature variation occurs along 
a line YjX = 2f/o/7o « 2 — const. This line crosses the critical line C = 0 at the 
roughening temperature 


, 2 , 1\[A 2 , 

k B T R = — 7 0 a 2 + U 0 = -Joa 

7T 7T 7T 


1 + 


he 

2 J 


(A9.48) 


with t c = 2\fAUol') 0 a i . Due to the potential Uo the stiffness 70 is modified to 
the effective one 7 e ir = 7 o(l + f c / 2 ), and correspondingly the microscopic coupling 
J = 7 0 a 2 /2 to the effective one J e a = %tra?/2. The roughening temperature is given 
as k B T R = 4J eff /ir. 

When the temperature is off from Tr., the parameter C is represented as 


4t(f c + 1) 


(A9.49) 


with t = (T - T r )/7r. For T < Tr, the potential strength Y initially decreases by 
renormalization until Y m = \JC, and then Y increases again. The final increment of 
y means the enhancement of the pinning potential U to fix the height to an integer 
value. The surface is smooth. The correlation length £ is determined as a length scale 
until y reaches to the order unity: 


, _ U XdY __ r°° 2 dY n 

fc “ ~ k 2{X - i)y ~ Jvc Jays/W^c ~ 2 


(A9.50) 


Thus the correlation length £ diverges at the roughening temperature T R as 


£ « exp \^{\t\t c ) 1/2 

r c 1 


~ exp 




(A9.51) 


The step free energy (3 is obtained as 


P ~ £ 1 ~ exp 


C 

VT^T. 


(A9.52) 


A9.6 Exact solution of surface roughening: Body-centered 
solid-on-solid (BCSOS) model 

There is an exact solution on the phase transition of the (100) surface of a body- 
centered cubic (bcc) crystal (Fig.A9.4a). There are eight nearest neighbors for each 
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« 



(d) 


Figure A9.4: (a) Body-centered cubic (bcc) lattice, (b) (100) face of the bcc lattice 
with the nearest neighbor bonds shown by dashed lines, and the 2nd nearest neighbor 
by solid lines, (c) Six possible configurations of surface heights, (d) The corresponding 
six vertices [187]. 
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atom in the bcc crystal. On looking down the (100) face nearest neighbor bonds 
form a square lattice (Fig.A9.4b). The solid-on-solid model assumes that the surface 
configuration is described by the height of the column on this square lattice. In a 
bcc crystal the nearest neighboring heights differ by odd numbers in a unit of half 
the lattice constant. When the nearest neighbor interaction marked by dashed lines 
in Fig.A9.4b is attractive and very strong, the height difference is restricted only 
to unity. Then for a square plaquette there are 6 possible configurations, as shown 
in Fig.A9.4c. When one draws a vertical line through the middle of an edge, and 
draw an arrow such that, when one follows the arrow direction, the right stays higher 
than the left. Then there are six vertex configurations as shown in Fig.A9.4d, and 
the model is equivalent to the six vertex model [129]. Among four edges from a 
vertex, two arrows flows inwards and the other two arrows go outwards. The model 
is originally introduced to explain the residual entropy of ice. In ice, each water 
molecule is surrounded by four neighboring water molecules by the hydrogen bonds. 
On the bond connecting two oxygen atoms, a hydrogen atom can be located near one 
of them. For each oxygen atoms, two hydrogen atoms lie close to it and two others 
stay away from it. The vertex in Fig.A9.4d represents the oxygen, and the arrow on 
the edge indicates that the hydrogen lies close to the oxygen. Thus the rule that the 
number of incoming and the outgoing arrows are the same is called the ice rule [153], 
We now impose a second neighbor interaction, —J x and —J v in x and y directions 
respectively, which are shown by solid lines in Fig.A9.4b. The interaction energies 
£i ~ ea for vertex configurations 1 to 6 are then given as 

£i = £2 = J y - J x , £ 3 = £4 = J z - J y , and £5 = £e = -J x - J r (A9.53) 

For an isotropic case, J x = J y — e/2 > 0, the configuration energies are e\ = £2 = 
£3 = = 0 and £5 = £e = —£. This is called the F model, and the exact solution is 

known [129]. The transition temperature is known to be 

= (In 2 )- 1 « 1.4427 • • • , (A9.54) 

and the singular part of the surface free energy behaves as 

F, ing * exp(-o|l - T/T r \- 1 ' 2 ), (A9.55) 

with a = 7 r 2 / 4 v / ln 4 . Step energies in [10] and [11] directions are also known as 

£[ 11 ] _ 2 y /2 J 1 

£ “ y/T^t\2 

OO 

+ ^(-l) n [l + (n - 1) tanh((n - 1)A) - n tanh(wA)] 

n=l 

£[ 10 ] _ 4 f 1 

x/r^r u 


£ 
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Figure A9.5: Step energies of BCSOS model in two orientations,[10] (solid line) and 
[11] (dotted line). For comparison, step energies of the Ising model is also shown. 
Dashed line in [10] direction and dash-dotted line in [11] direction. Crosses represent 
the difference of the step energies in two directions [187]. 

+ + (n- j)tanh((n - j)A) - (n - i)tanh((n - i)A)] j 

(A9.56) 

with 

A = -Inf+ 21n (l + \/l - f) (A9.57) 

and t — 4cxp(—2 e/hgT). At low temperatures, kgT < e, the step energy is 
anisotropic, £|u] > £[io|- At high temperatures, kgT > e, it becomes isotropic. Near 
the roughening temperature where t —> 1 and A « 2\/l — t « 2\/ln4 ^Tr/7 1 — 1, step 

energies behave as E|n] « Ep 0 ] « eA -3/2 exp(— 7r 2 /2 A) w exp(-a/^/l — T/T R ). They 
approach zero as T —♦ Tr — 0, as is shown in Fig.A9.5. The equilibrium shape of 
facets and the crystalline shape arc also obtained [92], 

Similar correspondence to the F model is found for the (100) and (110) faces of 
face centered cubic (fee) crystal [93]. 
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A14 Advance Rate of a Circular Step 


Detailed calculation of the advance rate of a circular step is given in the paper by 
Burton, Cabrera and Frank [43], and here I reproduce the calculation. 

The adatom density or the diffusion field around a circular nucleus with a radius 
p depends only on the radial variable r as u(r) = c(r)—c a 0 due to the symmetry, and 
the diffusion equation is written as 


d?u 1 du u 
dr 2 r dr x% 


(A14.1) 


in the stationary approximation. The solution which has no singularity at the origin 
and at infinity is written as 


u(r) = { 


v 


u(p) 


Ip(r/x s ) 

Iq{p/xb) 


u(p) 


Kg(r/x a ) 

Kg{p/x s ) 


for r < p 


for r > p . 


(A14.2) 


Here fo and K 0 are the modified Bessel function of the zeroth order. The i/-th order 
ones, I± v and K± v , satisfy the ordinary differential equation 


~d?w dw . o 
1 + +l,) ” = 0 ’ 


(A14.3) 


and I ± „(z) is finite at z —> 0 and K v {z) is finite at \z\ —> oo and |arg( 2 )| < |. For 
v = 0 , they can be expanded around z = 0 as 


h{z) = 
K 0 {z) = 





+ ••• 



+ 7 


h(z) + 


1 

(l !) 2 4 


+ < 1 + 5 > 


1 

w 


(A14.4) 



+ • • • (A14.5) 


with Euler’s constant 7 = 0.57721 56649 • • •. Their asymptotic forms for z 

I0(Z) ~ vfe[ 1 + ^ + '"l 

K o(z) ~ + 

By differentiation, the following relations hold; 


I^(z) = h(z), Kq(z) = -Ki(z), Io(z)K l (z) + I 1 (z)K 0 (z) = -. 


—» 00 are 
(A14.6) 
(A14.7) 


(A14.8) 
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The local equilibrium is expected around a circular nucleus; 

u(p) = c eq (p) - Coo. (A14.9) 

On the other hand, the advance velocity of the step is derived from the material 
conservation as 


v(p) 


, dc 
D a a 2 — 


dr 


p+ 


= D s a' 


Mp) 


h-oip/ %n) Iq(pI 3'h ) 


= -Dm 


2 «(P )_1 _ 

x s p/x„K 0 Io 


Ko(p/x„) I 0 {p/x H ) 
py i 
P Ko(p/x a )Io{p/x») 


= -Dm 


2 u (p) Ai/o + 7iA 0 
x 8 fC 0 Jo 

u(p). (A14.10) 


When the radius p is larger than x s , the approximation Ko(p/x a )Io(p/x s ) « x s /2p 
holds from Eq.(A14.6) and (A14.7), and n(p) is obtained as 


, , A, a 2 1 . . 2D a a 2 

v (p) = ; „ - /o , ( c °° _ c <*i(P)) = 


/?« 2 \ 
pksT ) 


p x s /2p ' . x a 

For the straight step with p —> oo, the advance velocity reduces to Vo as 
, V 2 D s a\ 

v(p -> 00) = - (Coo - C cq ) = VO. 

Xs 

Then the velocity of a circular step is approximated as 

v(p) = w o - ’ 


(A14.ll) 


(A14.12) 


(A14.13) 


where p c is the critical radius of the two-dimensional nucleation defined in (6.23) as 

Cc q (ia 2 _ ffa 2 


Pc = 


Coo ^eq fcflT Ap 


(A14.14) 
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A15 Advance Velocity of a Spiral Step 


There are various ways to represent the spiral curve. One way is to use an arc length 
s along the spiral and the angle 6 of the normal vector n from the y axis as shown in 
Fig.Al5.la. 

0 = 8(s). (A15.1) 


From Fig.A15.1a it is obvious that the radius of the curvature p is 
variation of the arc length ds and the associated angle change dff as 
the curvature is given as 

1 d0 


K = - 


P 


ds' 


related to the 
ds = pd6, and 

(A15.2) 


The same spiral can be represented by the polar coordinate (r, cj>) 
at the spiral center as 


4> = <P{r). 


with its origin 
(A15.3) 


Figure A15.1b shows that the variation of an arc length ds is written as 


ds = sj{dry + (rdcp) 2 = dr\J 1 + (rep 1 ) 2 , (A15.4) 


where the prime means the derivative by r: ' = d/dr. There is still another possibility 
of representing the spiral by the angle ip between the radial vector r and the tangential 
vector t of the spiral (Fig.A15.1c) [44]: 


ip = ip{r). 


(A15.5) 




(a) (b) (c) 

Figure A15.1: Three ways of representing a curved spiral, (a) By arc length s and the 
angle 6 of the normal vector n making with y axis, (b) The polar coordinate ( r,<p ). 
(c.) The length r of the radial vector r and the angle ip of the tangential vector t and 
r. Angles are related as ip = 6 + cp. 
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Figure A15.2: A spiral steadily rotating with an angular velocity u. v is the local 
velocity in the radial direction, and v is velocity normal to the spiral. 


As is clear from Fig.A15.1c, the three angles are related as 

ip = <p + 0. 

Also from Fig.A15.1e it is evident that the relation 

rip' = — tan ip 


holds. Then the curvature k is written as 

_ 1 _ <i0 __ dtp — dp 
P rfs dr^Jl + (rfi) 2 


sin ip 

-h ip cos ib. 

r 


(A15.6) 

(A15.7) 


(A15.8) 


When the spiral is advancing steadily, the spiral looks as if it is rotating steadily. 
By denoting the radial velocity at a point (r, <j>) as v, the point on the curve at a polar 
angle <p moves to 

r\ = r(t) + vdt (A15.9) 

after a time dt. Since the movement looks like a steady rotation with an angular 
velocity w, the angle <p(r lt t) on a spiral at time t should rotate to <p(r, t) after a time 
elapse dt: 

<p(ri,t) + u)dt = <p(r,t). (A15.10) 

(See Fig.Al5.2). Inserting (A15.9) into (A15.10) and expanding in terms of a small 
variable dt, radial velocity v is written as 


w 



(A15.ll) 
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The normal component of the step velocity should be 

v(r) = ucos ($ + 0 - ^ sin ip, (A15.12) 

and this velocity should coincide with the velocity given by the step curvature as 


v(r) = #o 1 - 


Pc 


(A15.13) 


From Eq.(A15.7-A15.8) and (A15.12-A15.13), the spiral ip = ip{r) should satisfy the 
first order differential equation 


dtp _ tttnip 1 
dr r p c cos ip 



(A15.14) 


Boundary conditions are set at r -> 0 and at r -* oc: At r —* oo, the spiral looks 
like a circle and the angle ip -* it j 2. At r — » 0, the radius of the curvature decreases, 
but it should be larger than p K , because if the radius is smaller than p c , the crystal 
melt back as is evident from Eq.(A15.13). Therefore, at r = 0 the curvature is p c . 
From (A15.8) ip should approach zero as r -> 0 in order to obtain a finite k. Then 
sin ip Jr ~ ip/r ~ dip/dr, and at r -» 0 the relation l/p c ~ 2dip/dr holds. Therefore, 
ip ~ r/2p c for small r. 

Since the angle ip satisfies the first order differential equation (A15.16) the solu¬ 
tion contains only one integral constant. Thus two boundary conditions cannot be 
satisfied in general. The solution exists only for some special values of Thus, 
Eq.(A15.14) is the eigenvalue equation to determine the eigenvalue w. From the 
numerical calculation, the eigenvalue is obtained to be 


u) « 0.33uo/pc. 


(A15.15) 


The period T for one turn of the spiral is T = 2ir/w. Since the advance velocity of 
the step is «o asymptotically at r —♦ oo, the step separation A = VqT is calculated as 

X = Hip « 19ft.. (A15.16) 

The step separation of a steadily growing spiral is 19 times larger than the critical 
radius p c of the two-dimensional nucleation. 

From the relation (A15.7), the polar angle <p for r -> 0 is determined from the 
relation r<p/ « -ip « -r/2p c to be 

(A15.17) 

2p c 

This is the Archimedes’s spiral turning right. The spiral turning left is written as 
<p = r/2p c , This is used to obtain the formula (15.1). 
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A21 Parabolic Coordinate 

The Cartesian coordinate (x, y, z) and the parabolic coordinate (£, i], (f>) are related 


as 


x - y/^n cos <p, y=^ s \n 0, z = ^(7?-£), 
where r = sjx 1 + V 1 4- z 2 = \{j] + ()■ Inverse relation is written as 

V 

£ = r + z, iy=r-z, tj>- arctan-. 

x 

The line element is represented in both coordinate as 

(■dsf = (dx) 2 + (dy) 2 + (dz) 2 = (h n dr]) 2 + (h{d() 2 + {h^d<j>f 

with 


= d 2 ^ 1 - *< 

The gradient vector of the field is 




(A21.1) 

(A21.2) 

(A21.3) 

(A2I.4) 


„ 5m_ dn „ a 
v * = 


4r; 3 m „ / 4£ 3u - jldu- 

7) + £ dr) 11 + V T] + (, + V d<t> 


(A21.5) 


with basis vectors (x,y,z) and (£,?), 0) in each coordinate system. The basis vector 
z is related to those of parabolic coordinate as 


1 

2 K 


■v- 


2 h, 


<■ 


(A21.6) 


The Laplacian is written as 


Am = 


p + C 


d ( dtt \ d ( du 
dv \9 t]) + 3£ \ 


+ 


1 d 2 u 


53- ( A2L7 ) 


d<t> 2 


By denoting the interface as p = r/,(£,<W, the unnormalized tangential vectors are 
written as 


*i = h n^fl + h ( £ 

. , % „ * 
t2 = h ’>W V + 


The normal vector is then obtained as 


n = IV 


J_. _ J_%x_ 1 dru ■ 

lO MC Kd<t> q 


= N 


4 rj 

V + Z 


n- 


4{ dr] 


i- 


(A21.8) 

(A21.9) 


1 dt] 
(A21.I0) 
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with N being determined from the normalization condition |n| = 1. 

The interface velocity of the dendrite in laboratory frame is written as v = Vz + 
rjh v ft. By using z and n given in Eq.(A21.6) and (A21.10), one obtains the interface 
normal velocity as 

(A2U1 > 

which is explicitly presented in Eq.(21.3). 
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A24 Dendritic Growth 


A24.1 Boundary integral equation for the diffusion-controlled 
growth 

In a steady state approximation, the diffusion equation in the moving frame is written 
as 

Lu = V 2 u + ~-^- = 0, (A24.1) 

*D C/Z f 

where z' is a new coordinate z' — z - vt and in = 2 Djv is the diffusion length. 
Hereafter in this section we drop the prime on 2 . By using the adjoint operator IJ 
defined by 

L\g(r^)~V\g-~^-, (A24.2) 

‘r> oz\ 

the Green’s theorem says that the integral over the region fli enclosed by the periphery 
Pi satisfies the relation: 


Idili [g(r,ri)Liii(ri) - tt(ri)L{(?(r,ri)] 


= idV, 


5(r ’ ri) ^r" tt(r,) ^fer + ^ ni ^ (r,ri)a(ri) 


(A24.3) 


When u(rG satisfies the diffusion equation (A24.1) and g is a Green’s function of the 
operator L\ as 

L{g(r,ri) = -<5(r - r,), (A24.4) 

then (A24.3) reduces to the integral equation over the crystal-liquid interface Psl- 


J dTifiig(T,Ti)~- = j dTi,si./»( r .ri)«(ri). 


(A24.5) 


Here we have used the boundary condition that far from the liquid-crystal interface, 
the diffusion field and its derivative vanish: u(r —> oo) = du/dn\ r -*<x> = 0. The 
Green’s function g is obtained by the Fourier transformation as 


g(r,rj) = J 


cpq 


-tq(r-ri) 


(2tt) 2 ql + (q z 




= — e 


f + l- n l 27T 


«-Ci )/ln K 


l r - r i| 

fa 


(A24.6) 


with K 0 being the modified Bessel function explained in Appendix A14. Another 
integral kernel h is obtained as 

h( r,rj) = ~ - y-n x , z g - c(r,)6(r - r,). (A24.7) 

on i to 


The coefficient c enters because the position r in Eq.(A24.5) lies on the boundary 
r SL . Instead of determining c(r) directly, one can impose a sum rule 
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on h as 

J dTi,SL*(r,ri) = -1. (A24.8) 

The condition is easily derived by inserting the trivial solution u = const of the 
diffusion equation in Eq.(A24.3) [30, 165]. Equation (A24.5) relates the diffusion 
field u to its normal derivative du/dn at the liquid-crystal interface. This relation is 
used in the numerical simulation described in Appendix A24.3. 

If there is a diffusion in the crystal with the diffusion constant D s , we can obtain 
another boundary integral equation similar to Eq.(A24.5) with the Green’s function 
g s and an integral kernel h$ . In them, the diffusion length is replaced by the diffusion 
length in the solid Id,s = 2£>s/v, and in hs the normal vector is replaced by n s = —n. 
Due to the difference of the direction of the normal vector, the sum rule for hs is 
different from (A24.8) and is written as [168] 

/dr liSL /i s (r,r 1 ) = 0. (A24.9) 

The diffusion field is continuous at the boundary; u s = «• When the heat is conducted 
in the crystal as well as in the liquid, the energy conservation boundary condition 
(18.7) is altered as 

= < A24 - 10 > 

In a symmetric model where D s = D, one gets simple relations; l DtS = Id, fls = g 
and h—hs = — <5(r—ri). By subtracting Eq.(A24.5) from the corresponding one in the 
crystal, and by utilizing the above relations and the energy conservation Eq.(A24.10), 
one gets 

/ar.^atr.r = l/ir liSl9 „ n 

= J dr 1>S L M r > r i) - h(r,ri)] u(ri) = «(r). (A24.ll) 

When the crystal grows steadily in z direction with the velocity vz, the normal velocity 
is v n = vn z and the integration can be transformed as / dl\si,n z = / dx\. Also by 
imposing the local equilibrium condition (18.11) or 

u(r) = A-da, (A24.12) 

one obtains a simple equation 

2 f°° 

A - dK = — / dx\g(r,ri) (A24.13) 

Id J—oo 

for the symmetric model. This is the relation often used in the theoretical analysis, 
as will be explained in Appendix A24.2. 
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A24.2 Microscopic Solvability Theory of Dendritic Growth 

The solvability theory is very mathematical and its detailed explanation is out of 
scope of the present notes. I sketch here briefly its main idea and the result [124, 105, 
33, 156], 

When the crystal is growing steadily, d(/dt = 0, and the diffusion constants 
in solid and melt is the same (symmetric model), the profile z = ((x) of a two- 
dimensional crystal is determined by (A24.13) or 

A — dn = ~ r dx i ))/,D Kq (. (A24.14) 
nln \ D J 

Here K 0 is the modified Bessel function of the zeroth order, explained in Appendix 
A14. Without capillarity (d = 0), Ivantsov parabola Civ(^) = —x 2 /2p satisfies 
Eq.(A24.14) at the supercooling A [154]. With the capillarity, the interface ((a;) 
deviates slightly from the Ivantsov parabola, Civ- By using the Ivantsov relation be¬ 
tween the supercooling A and the profile Civ and measuring the length in unit of the 
tip radius p of the Ivantsov parabola, one gets the relation 


<rA{x, ()k = I VP,*,C) - W.s.Gv). (A24.15) 


Here P = pjlp is the Peclet number, the function T 2 is defined by 


ra(P,ar,C) = \ /°° cte 1 e-' ,( «* ) -« I,)) A'o -1 ,) 2 + 

and a is the stability parameter defined by 

do __ 2 Dd 0 __ dov 
°~pP~ p 2 v ~ 2DP 2 ' 




(A24.16) 

(A24.17) 


The surface stiffness is assumed to have four-fold symmetry with the anisotropy factor 
A defined as 

A(x,0 = 1 — f cos 40 = 1 — e + (A24.18) 

e here means the strength of the anisotropy. Since a is proportional to the small 
capillary length do, it is also small. But the usual perturbation expansion in a is not 
applicable, since a couples with the highest derivative of the height or the curvature 


d 2 (,(x)/dx 2 

[1 + ( dC / ftr ) 2 ] 3 / 2 


(A24.19) 


One has to use singular perturbation method [16]. The situation is similar to the 
quantum mechanics [171]: In the Schrodinger equation 
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the Planck constant h = 2nh acts as a small parameter, but it is coupled to the highest 
second derivative of the wave function The approximation to the wave function 
'h can be obtained by a singular perturbation, known as WKB approximation [171]. 

Up to the linear order of the deviation of the interface profile from the Ivantsov 
parabola 

C(x) - Civ(x) = (1 + xy^Zix), (A24.21) 

Eq.(A24.15) is transformed to the following linear and inhomogeneous integro-differential 
equation [124] 


£Z 


a 


<PZ(x) (l + z 2 ) 1/2 
dx* + A(x) 1 ’ 

(1 + a; 2 ) 3 /' 1 r (x + a^i)(l + x\) 3 ! A 

2itA J 1 (m — mi)[l -f- mi) 2 ] 


Z(x 1 ) 


a 

(1 + a: 2 ) 3 / 4 ’ 


(A24.22) 


Here V denotes the principal value. For the existence of a nontrivial solution, the 
inhomogeneous term in the r.h.s. should be orthogonal to the zero eigenvector Z(x) 
of the adjoint operator £* as 

e <"-'>-/>(= < A2423 > 


This equation is an eigenvalue equation to determine the eigenvalue a. Z(x) is ob¬ 
tained by the WKB method, and 0(<r,e) is calculated to behave approximately as 


0(<r, c) ~ exp 



(A24.24) 


at small a and e. This equation shows that with an isotropic surface tension (e = 0), 
0 can never be 0 for finite a or finite velocity. To realize steady growth, an anisotropy 
e is necessary in the surface stiffness. For e / 0 there are infinitely many but discrete 
numbers of stationary solutions 

a & cr 0 c 7/4 (l + 2«) 2 ( with n = 0,1,2, * ■ •, 00 ). (A24.25) 


From the linear stability analysis around these stationary solutions, the steady state 
with the minimum a or the maximum velocity v is found to be stable, but the other 
solutions are unstable. 

The dendrite tip radius of the curvature p and the growth velocity v thus satisfy 
the scaling relation 

2 _ 2Ddo _ 2Ddo _7 ^ 


vp = 




a 


(A24.26) 
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as well as the Ivantsov relation (21.13), which reflects the total energy conservation 
and is influenced little by capillarity. The tip radius and the velocity is then deter¬ 
mined uniquely for a given supercooling A as 

P = ~ rfo-A-V 7 ^ (A24.27) 

ff 0 P{ A) o 0 

and 

u = A)^ 4 ~ ^AV' 4 . (A24.28) 

do TT^do 

In the last expression in Eq.(A24.27) and (A24.28) Peclct number -P(A) is approxi¬ 
mated by its limit for small supercooling in two dimensions. 

A24.3 Numerical Simulation of the Dendritic Profile 

Since the dendritic pattern is quite different from the flat interface, linear and weakly 
nonlinear analysis is insufficient. To find out the strongly nonlinear and far from 
nonequilibrium pattern, numerical simulation is a very useful method. There arc 
various different methods to simulated the profile of the crystal [135, 102, 103, 165, 
109, 110, 87]. The main problem of the simulation is the sharpness of the interface. 
If one solves the bulk diffusion equation numerically, a discrete mesh is necessary. 
When the interface moves, it will be off from the mesh points. Since the interface 
is susceptible to the fluctuation because of the Mullins-Sekerka instability, one has 
to be careful in tracing the interface motion. Also, from the microscopic solvability 
analysis, the anisotropy in the system is said to be essential to stabilize the dendrite 
tip. By using the mesh for the diffusion field, an artificial anisotropy from the mesh 
lattice can be introduced. There are various simulation methods to circumvent these 
obstacles. First we briefly summarize ideas of some proposed methods. 

By assuming the steady growth of the dendrite, the diffusion equation can be 
transformed to the integro-differential equation (A24.5) with the help of the Green’s 
theorem, as explained in Appendix A24.1. By solving this interface equation numer¬ 
ically there occur no problems stated above, since grid points lie on the interface. 
By integrating along one side of a steadily growing dendrite from the tail to the tip, 
one can search the selected velocity such as to produce a round tip which connects 
smoothly to the other side of the dendrite symmetrically [135, 102, 103]. This is the 
numerical realization of the solvability criterion of a symmetric dendrite in a station¬ 
ary code. Our simulation [165] relies on the integro-differential equation (A24.5), but 
it allows the variation of the system velocity so as to select the final profile sponta¬ 
neously. There arc other methods which don’t use boundary integral equation. In a 
phase field model, the interface is defined as the domain boundary of the phase field 
and is treated as a diffuse object, in order to get rid of the difficulty associated to the 
sharp interface [109, 110, 111, 117]. This is so far the sole model on which a three- 
dimensional simulation is performed [109, 111]. The phase-field model, however, uses 
a mesh to solve the diffusion equation, and the mesh-lattice anisotropy is expected 
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to affect the quantitative results. Recently a new algorithm with multiple meshes are 
proposed to solve the diffusion equation. These meshes are mutually shifted trans- 
lationally and rotationally [87]. With this multiple-mesh method the real dynamical 
evolution of the crystal interface is simulated in many systems. 

Here I describe our boundary integral algorithm applied first to study the selec¬ 
tion problem of a two-dimensional dendrite profile of a one-sided model, where the 
diffusion takes place only in the liquid phase [165]. The time evolution of the crystal- 
liquid interface is traced by solving the quasi-stationary distribution of the diffusion 
field by the boundary element method [30]. First assume that the profile TsL(t) and 
the frame velocity v(t) at time t are given. Then modification of Tsl and v to the 
asymptotic form is performed in the following processes. 

1. Calculate the diffusion field u at the boundary from the local equilibrium con¬ 
dition (A24.12), 

2. calculate the Green’s function g and the integral kernel h by using the frame 
velocity v , 

3. solve the linear equation (A24.5) for the normal gradient q = du/dn, since 
g, h , u and Fsl are known, 

4. the normal velocity is calculated by the conservation condition as v n (r) = 
~Dq{ r), 

5. then evolve the crystal profile in the normal direction n by r (t+dt) = r(t)+v n ndt 
to a new configuration Tsl(< + dt) at time t + dt, 

6. adjust the moving frame velocity v to the tip velocity of the dendrite n t i P by a 
rather ad hoc relaxation v(t + dt) = v(t) — (v(t) — v,; ]p (t))dt/r, 

7. and then iterate back to 1. and repeat the whole procedure until the steady 
state is realized. 

If the steady state is realized, the ad hoc relaxation of the frame velocity by the 
procedure 6. does not affect the resulting asymptotics. To solve the linear equation 
(A24.5) numerically by computer, one has to discretize the crystal front into a poly¬ 
gon, whose corner points are denoted by r ; - with j = 1, • • •, N. The diffusion field u 
and the normal gradient q = du/dn at a corner point r j are denoted by Uj and qj 
respectively. Every quantity at a point r on a polygon edge F, between r j and r ;+ i 
is interpolated linearly, for example, for u as 

n(r) = (0« j + <h(O u i+i (A24.29) 

and for the position r as 

r = 0,(Or J +02(Or J+ i (A24.30) 
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with an interpolation functions [30] 

MO = ~ MO = ^ (A24.31) 

with —1 < £ < 1. The integration along the edge I\, is written as 

ST dTj = I /-, ^ (A24.32) 

with Sj = ]r J+J — rj ( being the length of the edge segment iy Then the integral 
equation (A24.5) can be written as a matrix equation 

±G^=±H,u }f (A24.33) 

j=i j =l 


where 

Gij = j J t dig [Tj + y(0sj,r.] MO + TS I-i ^ ^ ^(^) s i-i* r «] MO 

(A24.34) 

and 

Hij = Ml fo + MOfy, r.'] MO + yp /, dih fo - MO^-i,^] MO 

(A24.35) 

with g and h given in Eqs.(A24.6) and (A24.7), respectively. The sum rule (A24.8) 
reduces to the condition 

= -1 (A24.36) 

;'=i 

and the diagonal element Hu is determined from 

= 1. (A24.37) 


The integration in Eq.(A24.34-A24.35) are performed by the four-point Gaussian 
quadrature method [1, 30], When the point r, happens to be one of the end points 
of the linear segment r,- where the integration is performed, the Green’s function g 
contains the logarithmic singularity as 


g(r,ri )« exp 



(A24.38) 


with Euler constant 7 . Integration then uses the Gaussian integration with a loga¬ 
rithmic singularity [1, 30]. 

Since our main interest lies in the evolution of the dendrite tip, we divide the 
interface into three parts along z: a tip region, a transition region and a tail region. 
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Evolution of the tip region is treated fully as described. The tail region consists of 
an Ivantsov parabola appropriate to the steady growth assumption with the frame 
velocity v. A transition region connects the tip and tail regions smoothly. 

Extension of the method to include the kinetic effect or dendrite tilting is straight¬ 
forward. When the kinetic coefficient K(n) is finite and anisotropic, the local equilib¬ 
rium assumption does not hold any more. Here one has to use the kinetic condition 
(18.8) [50]. In the steady state the diffusion field at the interface u(rj) is decomposed 
as 

U(ri) = A - dK -W)= Wo(ri) + W) Yn (A24 - 39) 

with 

«o(ri) = A - dn. (A24.40) 

The integral equation (A24.5) is now modified as 




D 


K{ n) 


dn\ = I < ^' l >SLhtt°( ri )- 


(A24.41) 


The kinetic effect introduces only a small modification in the integral kernel as g — 
Dh/K , which is readily calculated. The growth rate v„ of the interface is obtained 
from the normal gradient du/dn as in Eq.(18.7). Therefore, a simple modification of 
the integral kernel Gy can include the kinetic effect [50, 169]. 

In the directional solidification, the surface takes various periodic structures. In 
order to simulate an interface evolution with a periodic modulation with a periodicity 
A, one imposes a periodic boundary condition ((x + A) = ((x) in x direction. For 
numerical simulation, one period of interface is decomposed into N mesh points. 
Periodic images of the interface are taken into account in the kernel as 


E Gyq, = Z HyUj (A24.42) 

i =i i =i 

with 

oo oo 

Gy = E Gij+mN, and Hy = ff'j+m/v (A24.43) 

m=—oo m=—oo 

In practice, the influence of the ra-th image far from the system in consideration are 
negligible. We usually take the effect of images within the range of 51^ [167, 50]. For 
a dendrite in a channel with impermeable walls, the normal gradient of the diffusion 
field at the wall vanishes from the symmetry: du/dn = 0. This means the existence 
of mirror images at both walls: ((—x) — C(2A — x) — ((a') for a system within a 
channel of width A. Then the system has a periodicity 2A and one can appropriately 
modify Eq.(A24.42) in this case. 

In the unidirectional solidification, the temperature gradient is applied on the 
system. If orientations of the temperature gradient and of the crystalline axis are 
different, the tilted structure appears [2, 66]. In this case the tip of the dendrite 
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shifts transversally, and the locus of the tip is tilted from the temperature gradient 
by an angle <f>. The system is invariant in a moving frame with a transversal velocity 
v x = viiu\<p and a vertical pulling velocity v z = v. The diffusion equation in quasi¬ 
stationary approximation reads as 

2 (du du\ , , 

v " + s(s +, “ , 'sj =0 - (A24 - 44> 

The Green’s function g and an integral kernel h are now modified to 
<?( r , r i) = ^ ex P --— j- --- 

(A24.45) 



and 


h(r,ri) = 


1 (C _ Ci) ~ ( x ~ Xi)taxup 

2tt/d CXP -- 


-(n u + ni*tan 0 )A'o 


loCOS(p 


m • (ri - i 
|r — ri| cos 


iM£i£s)B <(, - tl) - < A24 - 46 > 


Further discretization and numerical simulation can be performed similarly as before 
[148]. 
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A26 Eutectic Growth Theory by Jackson and Hunt 


The surface supercooling AT at the eutectic crystal front is derived by Jackson and 
Hunt [91j. They considered both the lamellar and rod structure of an eutectic system, 
but here only the lamellar structure in two dimensions is summarized. 

The eutectic crystal is growing from the liquid with a concentration c^,. Con¬ 
centrations of the crystal a and p phases are cf and cf, the periodicity of lamellae 
is A, and the ratio of each phases are set r] and 1 — rj respectively. The material 
conservation requires the relation 

Coo = cg77 + cf(l - rj), (A26.1) 

and thus T} — (c§ ~ c x )/Ac and 1 - f) = (Coo - eg)/Ac with the miscibility gap 
Ac = Cg — c§ (> 0). 

The dimensionless concentration field u defined by 


u = 


c -c B 

Ac 


satisfies the diffusion equation in the liquid 


(A26.2) 


1 du _ 2 du 

ir c .dt~ U+ hdz 


= 0 . 


(A26.3) 


The last equality holds in the quasi-stationary approximation. The diffusion length 
in the liquid is defined by means of the chemical diffusion constant D c as 



(A26.4) 


In solid the material diffusion is assumed negligible (one-sided model). The thermal 
diffusivity is so high in all the phases that the temperature gradient is constant as 
Eq.(26.1): <?t = dT/dz = const or T(z ) = T E + Gy-z. 

The crystal-liquid interface is denoted a s z = ((x, t). Then the material conser¬ 
vation at the boundary yields 


du 

dz 


*=e 


where = (cf ,/S 
written as 


^ [u(:r, 0 - 


for 0 < x < r\\ in a phase 


(A26.5) 


2 

— [u(27, C) - «s] for < x < A in ft phase, 

— Ce)/Ac. The local equilibrium assumption at the interface is 


A T(x) = T e - T(C) = -GtC 

m a Acu(x, 0 + —k for 0 < x < 77A 

= « ° (A26.6) 

—mpAcu(x, C) + for 77A < x < A. 
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Here m„ j( j = \dTi,/dtf ' (t \ is the absolute slope of the liquidus lines with a and (3 crystal 
phases, k is the curvatures, 7 i, q and qr./j arc the surface tensions between the liquid 
and a ox ft phase, L a and Lp arc the latent heat of a and i3 phases respectively. 

From Eq.(A26.6) it is obvious that to obtain the supercooling AT at the average 
interface position ((), we have to know the average concentration field (u) and the 
average curvature (k). 

The field u is expected to be modified periodically with the same periodicity A of 
the interface modulation. Thus u can be expanded in a Fourier scries as 

oo 

u(x,z) = Uoc + £ B n c~ A ^~^ c-' 9 " 1 , (A26.7) 


where q n — 2irnjX. Since u has to satisfy the diffusion equation in the quasi-st&tionary 
approximation (A26.3), the decay rate A„ in z direction is found to be 

An = Zb 1 + \As 2 + (A26.8) 

The continuity equation (A26.5) is written at the average position {() as 


£A n B n e^ = ^ [ AU ( Z ) + £ B n e^ x ] , 


where 


A u{x) = | ^ 

Since the Fourier coefficient of Au(x) is calculated as 


-V 

-V 


for 0 < x < tjX 
for r]\ < x < A. 


/ dxAu(x)e 
Jo 


l/ x 

XJo 


Xq n 


V 2 


the coefficient B n for n ^ 0 is determined as 

4 ( ;-i-Q)W 2 s j ;l (A ?w „/2) 
lx>Xq n {A n — 2/lx)) 


B n = 


(A26.9) 

(A26.10) 

(A26.ll) 

(A26.12) 


Assuming the slow growth such that A <C In, the average concentration in front of 
the a or ft phase is calculated as 


1 f r i* 2 A 

(u n ) = — / u(j:,Odx- = M 00 + Bo + —P(i7) 

7} A JO 


2A 


<«0 = 


>>x 

(i 




2A 


where 


1 r* 

—r- / u(x, Qdx = Uoo + B 0 - — 

- n)x J r ,\ Zn(i -v) 

sm 2 (vijn) 


P(V), (A26.13) 


rv / X X ' 1 
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From Eq.(A26.13) one obtains 


(u a - vP) = 


c “ - ^ _ 2A P(rj) _ An P(r}) 


Ac 


or 


vAc — D t 


hv(l~r}) D c r)(l-rf) 

cl-4v( i-n) 


(A26.15) 


(A26.16) 


A/2 2 P(r)) ■ 

The growth velocity Eq.(A26.16) corresponds to the phenomenological expressions 
(26.3) and (26.5) with Q a in Eq.(26.5) at Coc = Cp being given as 


Q a = 


~ *l) 2 

2 P(v) 


(A26.17) 


In order to obtain the average curvature (k), we use the relation of the force 
balance at the triple point, where three phase, liquid, a and /?, meet. (See Fig.26.4). 
The surface tensions satisfy the relation (26.6a,26.6b) or: 


7lc» sin 6 a + 7 L/? sin Op = ^p 

7La cos 6 a - 7l^ cos Ofl = 0, 


(A26.18) 


where the angles 0„ and Op are defined in Fig.26.4. These contact angles give the 
average curvature as 


(«“) = \[ ndx =-^r f ~ds cos 0 --^r sin 6 a 
r/ A Jo pA J-e a ds pA 


j]\ 


V A 




(1 - p)A 


sin dp. (A26.19) 


Here we used the relation (A15.2) of the curvature k, the angle 6 of the normal vector 
and the arclength s: k = dO/ds. These curvatures are the same with those obtained 
in Eq.(26.7) by assuming that the interface is a part of circle. 

From Eq.(A26.6) the average supercooling at the interface, z = {() is expected 
to be constant for a and (3 phases as (A T a ) = = AT. Then the remaining 

unknown constant B 0 is determined as 


Moo + Bq — — - 


1 


m a + m 0 

+ 


2AT(p) (ma 
V 


Id 


nip 
1 -T) 


1 ( 27 LqT B 

AAc \ pL Q 


. a , 27lsTe . 

sin 0 a + —— - r sm Op 


(1 -v)Lp 


(A26.20) 


The interface supercooling AT is obtained as 


AT = 


m: 


—Ac 


2A T(p) 


+ m,g 1 l D p(l-p) 
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+ - 


(m- 1 + m l 

- s at - (£ + t) 




"1/3(1 - r i)Lp 


Here 


4 Ac P(??) y- 

AT m = Ti-1-n- 

*7(1 + m,p In 

is the minimum value of the interface undercooling, and 

K(1 - ??)sme„ + dpT) sin 9p] oc ~ 


(A26.21) 

(A26.22) 

(A26.23) 


is the corresponding periodicity of the lamellar structure. The capillary lengths d n 
and dp are defined as 


dn 


ll,gT R 

m a L„Ac 


and 


. _ l\.pTn 

13 nipLpAc' 


(A26.24) 


These equations (A26.21-A26.24) are more general than Eq.(26.10). 
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